{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyTorch Tests On/Off"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "from torch.autograd import Variable\n",
    "\n",
    "\n",
    "def _log_poisson_impl(n, lam):\n",
    "    normal = torch.distributions.Normal(lam, torch.sqrt(lam))\n",
    "    return normal.log_prob(n)\n",
    "\n",
    "\n",
    "class shapesys_constraint:\n",
    "    def __init__(self, nom_data, mod_data):\n",
    "        self.auxdata = []\n",
    "        self.bkg_over_db_squared = []\n",
    "        for b, deltab in zip(nom_data, mod_data):\n",
    "            bkg_over_bsq = b * b / deltab / deltab  # tau*b\n",
    "            self.bkg_over_db_squared.append(bkg_over_bsq)\n",
    "            self.auxdata.append(bkg_over_bsq)\n",
    "\n",
    "    def alphas(self, pars):\n",
    "        return np.product([pars, self.bkg_over_db_squared], axis=0)\n",
    "\n",
    "    def logpdf(self, a, alpha):\n",
    "        return _log_poisson_impl(a, alpha)\n",
    "\n",
    "    def expected_data(self, pars):\n",
    "        return self.alphas(pars)\n",
    "\n",
    "\n",
    "class Model:\n",
    "    def __init__(self, spec):\n",
    "        self.samples = {}\n",
    "        self.samples['sig'] = spec['sig']\n",
    "        self.samples['bkg'] = spec['bkg']\n",
    "        self.constraint = shapesys_constraint(spec['bkg'], spec['bkgerr'])\n",
    "\n",
    "    def expected_actualdata(self, pars):\n",
    "        return [\n",
    "            pars[0] * s + pars[1] * b\n",
    "            for s, b in zip(self.samples['sig'], self.samples['bkg'])\n",
    "        ]\n",
    "\n",
    "    def expected_auxdata(self, pars):\n",
    "        modparslice = slice(1, 2)\n",
    "        return self.constraint.expected_data(pars[modparslice])\n",
    "\n",
    "    def expected_data(self, pars, include_auxdata=True):\n",
    "        expected_actual = self.expected_actualdata(pars)\n",
    "\n",
    "        if not include_auxdata:\n",
    "            return expected_actual\n",
    "        expected_constraints = self.expected_auxdata(pars)\n",
    "        return np.concatenate([expected_actual, expected_constraints])\n",
    "\n",
    "    def logpdf(self, data, pars):\n",
    "        actual_data, auxdata = [data[0]], [data[1]]\n",
    "        lambdas_data = self.expected_actualdata(pars)\n",
    "\n",
    "        sum = 0\n",
    "        for d, lam in zip(actual_data, lambdas_data):\n",
    "            sum = sum + _log_poisson_impl(d, lam)\n",
    "\n",
    "        modauxslice = slice(0, 1)\n",
    "        thisauxdata = auxdata[slice(0, 1)]\n",
    "\n",
    "        modparslice = slice(1, 2)\n",
    "        modalphas = pdf.constraint.alphas(pars[modparslice])\n",
    "\n",
    "        for a, alpha in zip(thisauxdata, modalphas):\n",
    "            sum = sum + self.constraint.logpdf(a, alpha)\n",
    "        return sum\n",
    "\n",
    "\n",
    "def loglambdav(pars, data, pdf):\n",
    "    return -2 * pdf.logpdf(data, pars)\n",
    "\n",
    "\n",
    "def optimize(\n",
    "    objective,\n",
    "    init_pars,\n",
    "    par_bounds,\n",
    "    trf_to_unbounded=lambda x: x,\n",
    "    trf_from_unbounded=lambda x: x,\n",
    "    args=None,\n",
    "    constrain_index_value=None,\n",
    "):\n",
    "    data = args[0]\n",
    "\n",
    "    poi_index = None\n",
    "    constrained_var = None\n",
    "    if constrain_index_value is not None:\n",
    "        poi_index, poi_value = constrain_index_value\n",
    "\n",
    "        upto = init_pars[:poi_index]\n",
    "        fromon = init_pars[poi_index + 1 :]\n",
    "        nuis_init = np.concatenate([upto, fromon])\n",
    "\n",
    "        init_pars = nuis_init\n",
    "        constrained_var = Variable(torch.Tensor([poi_value]), requires_grad=True)\n",
    "\n",
    "    init_pars = trf_to_unbounded(init_pars)\n",
    "\n",
    "    #     print('init from: %s' % init_pars)\n",
    "\n",
    "    parameters = Variable(torch.Tensor(init_pars), requires_grad=True)\n",
    "    data = Variable(torch.Tensor(data))\n",
    "\n",
    "    optimizer = torch.optim.Adam([parameters])\n",
    "\n",
    "    def assemble_from_nuis(nuisance_pars):\n",
    "        tocat = []\n",
    "        if poi_index > 0:\n",
    "            tocat.append(nuisance_pars[:poi_index])\n",
    "        tocat.append(constrained_var)\n",
    "        if poi_index < len(parameters):\n",
    "            tocat.append(nuisance_pars[poi_index:])\n",
    "        return torch.cat(tocat)\n",
    "\n",
    "    for i in range(1000):\n",
    "        objpars = None\n",
    "        if constrain_index_value is None:\n",
    "            objpars = parameters\n",
    "        else:\n",
    "            objpars = assemble_from_nuis(parameters)\n",
    "        loss = objective(objpars, data, args[1])\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "    #         if i % 1000 == 0:\n",
    "    #             print(trf_from_unbounded(parameters).data.numpy())\n",
    "\n",
    "    result = parameters\n",
    "    if constrain_index_value is not None:\n",
    "        result = assemble_from_nuis(result)\n",
    "    result = trf_from_unbounded(result).data.numpy()\n",
    "    #     print('result %s' % result)\n",
    "    return result\n",
    "\n",
    "\n",
    "def unconstrained_bestfit(data, pdf, init_pars, par_bounds):\n",
    "    #     print 'fit', data, 'expected for init pars', np.array(pdf.expected_actualdata(init_pars))\n",
    "    result = optimize(\n",
    "        loglambdav,\n",
    "        init_pars,\n",
    "        par_bounds,\n",
    "        args=(\n",
    "            data,\n",
    "            pdf,\n",
    "        ),\n",
    "    )\n",
    "    return result\n",
    "\n",
    "\n",
    "def constrained_bestfit(constrained_mu, data, pdf, init_pars, par_bounds):\n",
    "    #     print 'fit', data, 'for pars', np.array(pdf.expected_actualdata(init_pars))\n",
    "    result = optimize(\n",
    "        loglambdav,\n",
    "        init_pars,\n",
    "        par_bounds,\n",
    "        args=(\n",
    "            data,\n",
    "            pdf,\n",
    "        ),\n",
    "        constrain_index_value=(0, constrained_mu),\n",
    "    )\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "spec = {\n",
    "    'obs': [55.0],\n",
    "    'sig': [10.0],\n",
    "    'bkg': [50.0],\n",
    "    'bkgerr': [7.0],\n",
    "}\n",
    "\n",
    "pdf = Model(spec)\n",
    "\n",
    "data = np.concatenate([spec['obs'], pdf.constraint.auxdata])\n",
    "init_pars = [1.0, 1.0]  # (mu, gamma)\n",
    "par_bounds = [[0, 10], [0, 10]]\n",
    "\n",
    "# constrained_bestfit(1.0, data, pdf, init_pars, par_bounds)\n",
    "# unconstrained_bestfit(data, pdf, init_pars, par_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "### The Test Statistic\n",
    "def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):\n",
    "    bestfit_nuisance_asimov = constrained_bestfit(\n",
    "        asimov_mu, data, pdf, init_pars, par_bounds\n",
    "    )\n",
    "    return pdf.expected_data(bestfit_nuisance_asimov)\n",
    "\n",
    "\n",
    "def qmu(mu, data, pdf, init_pars, par_bounds):\n",
    "    mubhathat = constrained_bestfit(mu, data, pdf, init_pars, par_bounds)\n",
    "    muhatbhat = unconstrained_bestfit(data, pdf, init_pars, par_bounds)\n",
    "\n",
    "    v_mbhh = Variable(torch.Tensor(mubhathat))\n",
    "    v_mhbh = Variable(torch.Tensor(muhatbhat))\n",
    "    v_data = Variable(torch.Tensor(data))\n",
    "    qmu = loglambdav(v_mbhh, v_data, pdf) - loglambdav(v_mhbh, v_data, pdf)\n",
    "    qmu = qmu.data.numpy()[0]\n",
    "    if muhatbhat[0] > mu:\n",
    "        return 0.0\n",
    "    if -1e-6 < qmu < 0:\n",
    "        log.warning('WARNING: qmu negative: %s', qmu)\n",
    "        return 0.0\n",
    "    return qmu\n",
    "\n",
    "\n",
    "def pvals_from_teststat(sqrtqmu_v, sqrtqmuA_v):\n",
    "    from scipy.stats import norm\n",
    "\n",
    "    CLsb = 1 - norm.cdf(sqrtqmu_v)\n",
    "    CLb = norm.cdf(sqrtqmuA_v - sqrtqmu_v)\n",
    "    CLs = CLb / CLsb\n",
    "    return CLsb, CLb, CLs\n",
    "\n",
    "\n",
    "def runOnePoint(muTest, data, pdf, init_pars, par_bounds):\n",
    "    asimov_mu = 0.0\n",
    "    asimov_data = generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds)\n",
    "\n",
    "    qmu_v = qmu(muTest, data, pdf, init_pars, par_bounds)\n",
    "    qmuA_v = qmu(muTest, asimov_data, pdf, init_pars, par_bounds)\n",
    "\n",
    "    sqrtqmu_v = np.sqrt(qmu_v)\n",
    "    sqrtqmuA_v = np.sqrt(qmuA_v)\n",
    "\n",
    "    sigma = muTest / sqrtqmuA_v if sqrtqmuA_v > 0 else None\n",
    "\n",
    "    CLsb, CLb, CLs = pvals_from_teststat(sqrtqmu_v, sqrtqmuA_v)\n",
    "\n",
    "    CLs_exp = []\n",
    "    for nsigma in [-2, -1, 0, 1, 2]:\n",
    "        sqrtqmu_v_sigma = sqrtqmuA_v - nsigma\n",
    "        CLs_exp.append(pvals_from_teststat(sqrtqmu_v_sigma, sqrtqmuA_v)[-1])\n",
    "    return qmu_v, qmuA_v, sigma, CLsb, CLb, CLs, CLs_exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XlU1Nf5+PH3nRkY9kU2WUXAHRSUAYwKLkRB4opEE7ekidakSZMmX5ukbdI2vyZpkia1aWsSY02r0aiIiHFDoxijMe5GDcFdFNxFcUMR+fz+wCGoLDPDIDPDfZ3Tc4T5zOc+k3P6zOXe+3keoSgKkiRJkm1RNXcAkiRJkvnJ5C5JkmSDZHKXJEmyQTK5S5Ik2SCZ3CVJkmyQTO6SJEk2qMHkLoSYJYQ4K4TYV8frQgjxkRDikBBijxCiu/nDlCRJkoxhyMz9v0BKPa+nAu3u/G8y8HHjw5IkSZIao8HkrijKBqCknkuGAbOVKt8DHkIIf3MFKEmSJBlPY4Z7BAInavxcdOd3p+69UAgxmarZPc7Ozj06duxo9GAnT56s/rdGo8He3h57e3u0Wi1qtdro+0mSJFmTHTt2nFcUxaeh68yR3A2mKMoMYAZAbGyssn37dqPvsWHDBvLy8gDw8fHh6tWrlJWVAeDm5kZoaCht2rShTZs2tGrVCiGE+T6AJElSMxNCFBpynTmSezEQXOPnoDu/axIJCQnk5eWhKArnzp1j0KBBhIWFcezYMQoLCzl8+DB79uwBwMXFpTrZh4aG4uXlJZO9JEktgjmS+1LgOSHEfCAeKFUU5b4lGXPRL8HcvHmTCxcukJuby+DBg4mLiyMuLg5FUbhw4UJ1sj927Bj79lUd9HF2dq6e1YeGhuLj4yOTvSRJNqnB5C6E+BLoC3gLIYqAPwJ2AIqifAKsAAYDh4DrwJNNFaxeQEAAR48excHBgYqKClasWIFaraZ79+4IIfD29sbb25vY2FgURaGkpITCwsLqZJ+fnw+Ak5NTdbJv06YNfn5+MtlLkmQTGkzuiqI81sDrCvArs0VkgA4dOnD06FGcnZ3JyclhypQpfPXVV6jVarp163bXtUIIvLy88PLyonv37iiKwqVLl+5K9j/99BMADg4OdyX71q1bo1LJ57wkSbI+D3RD1VzCwsIAUKlU6HQ6Vq5cyciRI8nJyUGlUhEVFVXne4UQeHp64unpSXR0NAClpaXVyziFhYXs378fAK1WS0hISPUyjr+/v0z2kiRZBatM7t7e3qjVam7fvk1AQADLly9nypQpVFZWkp2djVqtpnPnzgbfz93dnW7dulXP+i9fvlyd6AsLCzl48CBQtd4fFxdHUlISGo1V/qeTJKmFsMoMJYTA19eXU6dOodFoGDBgAC+//DJbt25l/vz5ZGVloVKpMOUcPVQdqYyKiqr+C+Dq1asUFhZSUFDAxo0bKSgoYNiwYQQFBZnzY0mSJJmN1a4xtG3bFgBPT0/CwsI4fvw4//3vfxk7diz+/v5kZmZWz7gby8XFhS5dupCens7YsWMpLy9n1qxZrF69mlu3bpllDEmSJHOy2uQeHFx1tP7q1asIIRg1ahRvvPEG169fZ9y4cfj5+bFgwQIOHz5s1nEjIiJ49tln6d69O5s3b+bTTz/l+PHjZh1DkiSpsaw2ueuXRMrKyggNDaVdu3ZcuXKFP//5zzg4ODB+/Hi8vb2ZP38+R48eNevYWq2WRx55hPHjx3P79m0+//xzVq1aJWfxkiRZDKtN7i4uLri6ugJVG6zl5eU8++yz/Pvf/6agoABHR0fGjx+Pp6cnX375JYWFBj2xa5SwsDCeeeYZdDodW7Zs4eOPP+bYsWNmH0eSJMlYVpvcoWppRqVSUVJSgp+fH2FhYTg5OfHSSy8BVU+kTpgwAXd3d+bNm8eJEycauKPx7O3tGTx4MBMnTgTgf//7HytWrKC8vNzsY0mSJBnKqpN7YGAglZWVFBYWotPpKCkp4fe//z0rV65kxYoVQNUMf8KECbi4uDB37ty7qkqaU2hoKFOmTCE+Pp5t27bx8ccfc+TIkSYZS5IkqSFWn9wBKioqcHV1xcnJiaCgINq3b89LL71UvQbu6urKhAkTcHR0ZM6cOZw61TSlb+zt7UlJSeHJJ59EpVIxZ84cli1bxs2bN5tkPEmSpLpYdXIPCAgAqs69Hzt2jNjYWA4ePMjbb7/N/v37+fe//119rbu7OxMnTkSr1TJnzhzOnDnTZHGFhIQwZcoUevbsyc6dO5k+fbrZT+1IkiTVx6qTu52dHa1bt0ar1XL48GFiY2NRqVQ4OzszcOBA/vznP3P+/Pnq6z08PJgwYQIajYbZs2dz7ty5Jo1t4MCB/OIXv8De3p4vvviCnJwcbty40WRjSpIk6Vl1coeq2futW7c4e/YsiqLQpUsXdu3axXvvvceVK1d444037rq+VatWTJw4EZVKxezZs7lw4UKTxhcUFMQvf/lLevXqxQ8//MD06dM5cOBAk44pSZJk9ck9KCiI27dvA3D48GESEhIoLy/nxo0bPPPMM3z66afs3bv3rvd4eXkxYcIEKisr+d///kdJSX0tYhtPo9GQnJzM008/jaOjI19++SVLliyp7iAlSZJkblaf3PWbqlqtliNHjhAQEEBwcDBbtmzhjTfewN3dnRdffJGqysQ/8/HxYcKECVRUVDB79mwuXbrU5LEGBAQwadIk+vTpw549e5g+fToFBQVNPq4kSS2P1Sd3b29v7O3tcXFx4fDhw1RWVhIfH8/FixcpKSnhzTffZN26deTk5Nz3Xj8/P8aPH8/Nmzf53//+R2lpaZPHq9Fo6N+/P5MmTcLZ2ZkFCxaQlZXF9evXm3xsSZJaDqtP7iqVioCAAG7fvk1ZWRmnTp2iU6dOuLm5sWXLFqZMmULnzp35v//7v1qPJPr7+zNu3DjKysqYPXs2V65ceSBx+/v7M2nSJPr27Ut+fj7Tp0+v7hAlSZLUWFaf3KFqaUY/6z58+HB1E4+jR49y4cIF/v73v3P48GH+8Y9/1Pn+sWPHcvXqVWbPns3Vq1cfSNxqtZqkpCQmT56Mm5sbmZmZZGZmcu3atQcyviRJtssmkntQUBCKouDl5VV9nrxHjx5oNBq2bNnCwIEDeeSRR/jLX/5S5/n24OBgHn/8cUpLS5k9e/YDTbB+fn489dRT9O/fn/379zN9+nT27dt33z6BJEmSoWwiues3Vd3c3Dhx4gQ3btzA0dGRbt26sXfvXq5fv84HH3zAjRs3+P3vf1/nfdq0acNjjz3GxYsXmTNnzgM9zaJWq+nTpw+TJ0/G09OTrKwsFi5c+MD+ipAkybbYRHJ3dXXFzc0NAEVRqkv8xsfHU1FRwY4dO2jfvj3PP/88s2bNYteuXXXeq23btowZM4bz588zZ86cB/7Qka+vL7/4xS9ITk7m4MGDTJ8+nT179shZvCRJRrGJ5A5Vs/eSkhLs7e2rl2Z8fHwICwtj27Zt3L59m9dffx0vLy9eeOGFepNleHg4o0eP5syZM3zxxRcPvDaMSqWiV69eTJkyBS8vL7Kzs8nMzKSiouKBxiFJkvWyqeReWlpKcHAwhw8frk7eCQkJXLlyhZ9++gkPDw/eeustvv32WzIzM+u9X7t27cjIyODUqVPMnTu3WUr4ent78+STT5KcnMxPP/1EVlZW9QNbkiRJ9bGZ5K7vzOTp6cmlS5eqnzqNiIjAy8uL77//HoCnnnqKbt26MXXq1AbX1Dt27Eh6ejpFRUXMmzevWTot6WfxqampFBQUsGTJEiorKx94HJIkWRebSe7+/v4IIRBCAHDo0CGgqmJkXFwcxcXFFBUVoVarmTZtGsePH+eDDz5o8L6dO3dmxIgRHD9+nPnz5zdbK724uDgGDBjAvn37WLZsmVyDlySpXjaT3O3t7fH19aWkpARPT8+7Sux269YNrVbLli1bAOjbty8jR47knXfeobi4uMF7R0VFMWzYMI4cOcKaNWua7DM0pHfv3iQmJrJr1y5WrVolE7wkSXWymeQOVevuxcXFhIeHc+zYseoNSK1WS0xMDPn5+Vy+fBmA999/n4qKCl577TWD7t2tWzd0Oh3bt29v0lrwDenbty8JCQls3bqVdevWNVsckiRZNptL7jdu3MDX15dbt27d1TM1Li4ORVHYtm0bUNXc+qWXXmLOnDnVM/qG9OvXDwcHh2adNQshGDhwID169GDjxo1s2LChWeKQJMmy2VRy12+qCiFQqVTV6+5QtdHaoUMHduzYUb1u/rvf/Y7WrVvzwgsvGLRJ6ejoSP/+/Tl27Fiz1oERQpCWlkbXrl3Jy8ur3iyWJEnSs6nkrq8Qefbs2eojkTXFx8dTVlZWXd/d1dWVd955hy1btjBv3jyDxujevTutW7dm9erVzba5ClUJftiwYXTu3Jnc3Fx27NjRbLFIkmR5bCq56ytE6tfdz5w5c9fj+23atMHPz48tW7ZUL6tMmDCB2NhYXn31VYPqyahUKlJTU7l8+TIbN25sss9iCJVKxciRI2nXrh3Lli1jz549zRqPJEmWw6aSO1Stu58+fZrQ0FCAu2bvQgji4+M5e/Ysx44dA6oS5LRp0yguLubdd981aIyQkBCioqLYtGkTFy9eNPdHMIparebRRx+lbdu2LFmyRJYNliQJsNHkrl8/d3Jyum9pJioqCicnp7s2UXv16sWYMWN4//33KSwsNGic5ORkVCoVq1evNl/wJtJoNIwZM4agoCCysrI4ePBgc4ckSVIzs7nkrt9U1S/N1CxFAFWJMDY2lv3799/VO/Xdd99FCMErr7xi0Dhubm706dOHgoKC+75AmoO9vT2PP/44fn5+LFiwoLp4miRJLZPNJXdXV1dcXV2rk/v169c5derUXdfExsaiUqnYunVr9e9CQkKYOnUqCxYsMHgtvWfPnnh6erJq1SqLqPni4ODAuHHj8PLy4ssvv7zrKKgkSS2LQcldCJEihNgvhDgkhHi1ltdDhBB5QohdQog9QojB5g/VcEFBQdXJHbhvZu3q6kpkZCS7du26q+Ljb3/7WwIDAw0+GqnRaBg0aBDnz5+vPj/f3JycnBg/fjxubm7MnTuXkydPNndIkiQ1gwaTuxBCDfwbSAU6A48JITrfc9kfgIWKosQAY4Dp5g7UGIGBgVy8eBEhBK1bt6512SQ+Pp7y8vK7ars7Ozvz3nvvsXPnTv773/8aNFb79u2JiIhg/fr1FtMez8XFhfHjx+Pg4MAXX3zB2bNnmzskSZIeMENm7nHAIUVRjiiKUg7MB4bdc40CuN35tzvQrNNFfWemkydPEh4ezokTJ+6ryR4QEEBwcDBbt269a5b+2GOP0bNnT373u99VlyqojxCCQYMGcevWLdauXWveD9II7u7uTJw4EY1Gw+zZs7lw4UJzhyRJ0gNkSHIPBGou3hbd+V1NfwLGCSGKgBXA87XdSAgxWQixXQix/dy5cyaEa5iAgACEEBQVFREeHk5lZWX10cea4uPjuXjx4l2nS4QQ/OMf/+DMmTO8/fbbBo3n7e1NQkICu3btMqgQ2YPi6enJhAkTUBSF2bNnc+nSpeYOSZKkB8RcG6qPAf9VFCUIGAzMEULcd29FUWYoihKrKEqsj4+PmYa+n729PT4+PhQXFxMcHIydnd1dpQj0OnXqhJub2321ZXQ6HRMmTODvf/+7wSdhEhMTcXFxsbhqjd7e3owfP57y8nJmz55t0F8jkiRZP0OSezEQXOPnoDu/q+kpYCGAoiibAQfA2xwBmkpfIVKtVhMaGlprklapVOh0Oo4ePXpfpcd33nkHOzs7pk6datB4Wq2W5ORkioqK+OGHH8zyGcyldevWjBs3jmvXrjFnzhyL2RuQJKnpGJLctwHthBBthRD2VG2YLr3nmuPAAAAhRCeqknvTrbsYICgoiBs3blBSUkJ4eDgXL16861y7Xo8ePdBoNPfN3gMCAnjttdfIzs42uLRu165dCQoK4uuvv37gfVcbEhgYyOOPP86lS5eYM2dOg12oJEmybg0md0VRKoDngFzgJ6pOxfwohHhTCDH0zmUvA5OEED8AXwJPKM28NqHfVC0uLiYiIgK4/0gkVFV67NatG3v37uX69et3vfbSSy/Rpk0bXnzxRYOaUwshSElJ4dq1a3zzzTdm+BTm1aZNG8aMGcP58+eZO3euxX0BSZJkPgatuSuKskJRlPaKooQrivLWnd+9oSjK0jv/zlcUpZeiKN0URYlWFKXZn8n38fHBzs6OoqIiWrVqhYeHR53r5/Hx8VRUVNxXWdHR0ZG//e1v7N27l5kzZxo0bmBgIDExMWzZsoXz5883+nOYW3h4eHXj7+bqCytJUtOzuSdU9WpWiBRCEB4eztGjR2t9ktTHx4ewsDC2bdt23+vp6ekkJiby+uuvG3zaZMCAAdjZ2Vnc5qpehw4dGDFiBCdOnGDBggUG/VUiSZJ1sdnkDj9XiKyoqCA8PJzy8vI6H8lPSEjgypUr91VVFEIwbdo0Lly4wJtvvmnQuM7OzvTt25fDhw9z4MCBRn+OphAZGcnQoUM5fPgwixYtsojyCZIkmY9NJ/egoCAqKys5ffo0bdu2RQhR59JMREQEXl5etbbci4mJ4amnnuKf//wn+/fvN2hsnU6Hj48Pubm5Fjszjo6OZvDgwezfv5/s7GyDSi5IkmQdbDq519xUdXBwqLU7k54Qgri4OIqLiykqKrrv9b/85S84ODjw5z//2aCx1Wo1KSkpXLx4kc2bN5v+IZqYTqcjOTmZH3/8ka+++soil5EkSTKeTSd3Nze36gqRULWZeOrUqTrPeUdHR6PVamudvfv5+TFlyhQWLFhg8INNYWFhdOrUiW+//daiHx7q1asXSUlJ7N69m5UrV8oEL0k2wKaTO1TN3vUz8bqqROrZ29sTExNDfn5+rcn4N7/5DRqNhr/97W8Gjz9w4EAURWHNmjUmRP/gJCUl0bNnT7Zt28bXX38tE7wkWbkWkdwvXrzI9evX8ff3x9HRsd6Zd1xcHIqi1FrCNyAggCeeeILPP/+c06dPGzS+h4cHDz30EPv27TO4y1NzEELw8MMPExsby3fffceGDRuaOyRJkhrB5pN7zc5MKpWq1u5MNXl6etKhQwd27NhR6xnwqVOncuvWLf7+978bHEPv3r1xc3Nj5cqVFr1pKYRg8ODBREdHs379er777rvmDkmSJBPZfHL39/cHuGvd/dq1a/fVkqkpPj6esrIy9u7de99rERERZGRk8PHHHxt87t3Ozo6BAwdy5swZdu7cacKneHCEEAwZMoQuXbqwZs0ai2lCIkmScWw+uWu1Wnx9fe9K7kCtVSL12rRpg5+fH1u2bKl1hv/qq69y5coVpk83vCdJ586dCQ0NZd26dRZf10WlUjFixAjat2/PihUr7jv7L0mS5bP55A4/V4hUFAVXV1d8fX3rXXcXQpCQkMDZs2drrQMfHR1Namoq06ZNu68eTX33TElJ4caNG+Tl5Zn6UR4YtVpNRkYGgYGB5OTk1Fp0TZIky9VikntZWRkXL14Eqmbvx48fp7y8vM73REZG4uTkxPfff1/r66+++irnzp1j1qxZBsfh5+eHTqdj+/btBm/INieNRkNGRgZqtZqFCxfKOjSSZEVaTHIHqo9ERkRE1NmdSU+j0RAbG8uBAwdqnbX26dOHhx56iPfff9+opNe3b18cHBwstu7Mvdzd3RkxYgRnzpxh5cqVzR2OJEkGahHJ3dfXFzs7u+p195CQEDQaTYMPI8XGxqJSqdi6det9rwkheO211zh+/Djz5883OBZHR0cGDBhAYWEhP/74o3EfpJm0a9eO3r17s2vXLotrRCJJUu1aRHKvWSESqmbldXVnqsnV1ZXIyEh27dpVa+3ztLQ0IiMj+etf/2rUEceYmBhat27NmjVr6l0asiT9+vWjTZs2LF++nLNnzzZ3OJIkNaBFJHe4u0IkVK27X7hwocHjjPHx8ZSXl7Nr1677XhNC8Oqrr5Kfn89XX31lcCwqlYrU1FQuX77Mxo0bjfsgzUSlUpGeno69vT2ZmZlW86UkSS1Vi0rut2/frj7fru/OVN+RSKh6KjU4OJitW7fWOjsfPXo0bdu25Z133jFqDT0kJISoqCi+++676o1eS+fq6kp6ejoXLlxg2bJlVrFnIEktVYtK7vDzpqqXlxfu7u4GFQGLj4/n4sWLHDx48L7XNBoNU6dOZcuWLUa31ktOTkalUrF6dbM3rjJY27Zt6du3L3v37r2vc5UkSZajxSR3Nzc3XFxcqtfdhRCEhYXV2Z2ppk6dOuHm5lZrtUiAJ598Ej8/P9555x2jY0pMTKSgoMDgSpOWoE+fPkRERLBq1SpOnjzZ3OFIklSLFpPchRAEBQVVJ3eoWpq5efPmXb+rjUqlQqfTcfTo0VrLFjg4OPDiiy+yevVqo2ezCQkJtGrVilWrVllNNyQhBCNGjMDZ2ZnMzExu3LjR3CFJknSPFpPcoWpppqSkpPrxf313pobW3QF69OiBRqOp9VgkwDPPPIObmxt//etfjYpJo9EwaNAgzp8/X+e9LZGTkxOjRo3i8uXL5OTkyPV3SbIwLS65w89FxBwdHQkMDDRoScTR0ZGoqCj27t1b60zV3d2dX/3qV2RlZRndN7V9+/a0a9eO9evXc/XqVaPe25yCg4N5+OGHKSgosOhuU5LUErWo5B4QEABwVxu98PBwTp48aVCNGJ1Ox61bt9i9e3etr7/wwgtotVree+89o2MbNGgQFRUVrF271uj3Nqf4+Hg6derE119/zfHjx5s7HEmS7mhRyV2r1eLj43PXJqD+SOSRI0cafL+/vz9BQUFs37691mUIPz8/fvGLXzB79uxa+7DWx8vLi4SEBHbv3t3gHoAlEUIwdOhQPDw8WLRoUZ0tDCVJerBaVHKHn9vu6ZNzQEAADg4OBp9W0el0XLhwgaNHj9b6+tSpU6msrOTDDz80OrbExERcXFysro+pg4MDGRkZXL9+nezsbKuKXZJsVYtM7jUrRKpUKsLCwurtzlRT586dcXJyqrOJRWhoKI899hgzZszgwoULRsWm1WpJTk6muLjY6mq4+Pv7k5qayuHDh2WLPkmyAC0uuddsu6cXHh7OlStXDKqZotFoiImJYf/+/ZSWltZ6zSuvvMK1a9f417/+ZXR8Xbt2JSgoiK+//trqjhh2796drl27sn79eoOWuSRJajotLrnfWyESfl53N3RpJjY2FkVR6jzTHhkZyZAhQ/joo4+MPv0ihCA1NZVr164Z/cRrcxNCkJaWhre3N4sXL+bKlSvNHZIktVgtLrmrVCr8/f3vSu5ubm74+PgYnNw9PDxo3749O3furC5Edq/XXnuNkpISPvvsM6NjDAgIoHv37mzdupVz584Z/f7mZG9vz6OPPkp5eTlZWVkW3RBckmxZi0vuULXufurUqbueCA0PD6ewsNDgxhs6nY5r167x008/1fp6z549SUpK4oMPPjCpgmL//v2xs7MjNzfX6jYofXx8eOSRRygsLGTdunXNHY4ktUgtNrnfvn37rlZ34eHh3L59u97uTDWFh4fTqlWrOjdWoWr2XlxczBdffGF0jM7OzvTr14/Dhw+zf/9+o9/f3Lp27Ur37t3ZtGmT0Q91SZLUeC0yude2qdqmTRuDujPpCSGIjY3lxIkTdfZDHThwIDExMbz77rsm1Y2JjY3Fx8eH3NzcOpd/LFlqaiqtW7cmOzu7wbr5kiSZV4tM7vdWiASws7OjTZs2RlVnjI6ORqPR1Dl71zfzOHDgANnZ2UbHqVarSUlJ4dKlS3z33XdGv7+56RtsK4pCZmamVX5BSZK1apHJXQhBYGDgfU+ChoeHc/78+TqPON6roXozAOnp6bRr146//vWvJq2dh4WF0alTJ7799lurnP22atWKYcOGcfLkSdasWdPc4UhSi2FQchdCpAgh9gshDgkhXq3jmkeFEPlCiB+FEPPMG6b5BQYGcuHCheoKkVCV3MHwI5Hwc72Zuh46UqvV/Pa3v2XHjh18/fXXJsU6aNAghBCsXLnSpPc3t06dOpGQkMDWrVutpim4JFm7BpO7EEIN/BtIBToDjwkhOt9zTTvgNaCXoihdgBebIFazurdCJFSd8nB1dTUquevrzWzbtq3Omfn48eMJCAgwupmHnru7O0lJSRw4cICCggKT7tHckpOTCQoKYunSpUY/uStJkvEMmbnHAYcURTmiKEo5MB8Yds81k4B/K4pyEUBRlIYf9WxmtSV3IQTh4eEcOXLEqPPZsbGx9dab0Wq1vPTSS+Tl5dXZzakhCQkJ+Pr6smrVKqtsTq1Wqxk1ahRqtZqFCxcafORUkiTTGJLcA4ETNX4uuvO7mtoD7YUQm4QQ3wshUmq7kRBishBiuxBie3M/nKOvEHnvuntERAQ3btwwqjJjly5d6q03AzB58mQ8PT2Nbuahp1arSUtLo7S01Gprt7i7uzNy5EjOnj3LihUrmjscSbJp5tpQ1QDtgL7AY8BnQgiPey9SFGWGoiixiqLE+vj4mGlo0+k3VWsup4SFhSGEMGppxpB6M66urjz//PMsWbKE/Px8k+INCQkhOjqazZs3G1QHxxJFRETQp08fdu/eXWddfEmSGs+Q5F4MBNf4OejO72oqApYqinJLUZSjwAGqkr1FCwwM5Pr163edQnF0dCQgIMDohtUN1ZsBeP7553FycuLdd981OeaHH34YrVbLihUrrO7JVb2+ffsSGhrK8uXLa+1JK0lS4xmS3LcB7YQQbYUQ9sAYYOk91yyhataOEMKbqmUaiy8LWNu6O1SdmikuLr7rJE1DatabqeuBJW9vbyZNmsS8efMoLCw0KWYnJyeSk5MpLCy0urLAeiqVivT0dBwcHMjMzOTmzZvNHZIk2ZwGk7uiKBXAc0Au8BOwUFGUH4UQbwohht65LBe4IITIB/KAqYqiWPyRCD8/PzQazX1dk8LDw1EUxeiytfp6M/Utu7z88ssIIfjggw9MihkgJiaGoKAg1qxZY9QXkCVxcXEhPT2dkpISli1bZrV/hUiSpTJozV1RlBWKorRXFCVcUZS37vzuDUVRlt75t6IoykuKonRWFCVKUZT5TRkRv2G6AAAgAElEQVS0uahUKgICAu6buQcFBaHVao1emtHXm9m+fXud1wQHBzNu3DhmzpxpcsVHIQSPPPIIZWVlJp+dtwShoaH069ePffv21fvfTJIk47XIJ1RrCggIuK9CpLHdmfT09WaOHz9e71ryb3/7W27cuME//vEPk+P28/MjPj6enTt3Gt2v1ZL07t2bdu3akZube1dvW0mSGqfFJ/egoCBu3759XzIODw/n8uXLnD9/3qj76evNbN26tc5rOnbsyIgRI/jXv/7F5cuXTYobqjYmXV1dWbZsmdXWTRdCMHz4cFxcXMjMzLTaZSZJsjQtPrnXt6kKcOjQIaPu5+joSGRkZL31ZqCqHHBpaSmffvqpkRH/TKvVkpKSwpkzZ+r9MrF0Tk5OjBo1isuXL5OTkyPX3yXJDFp8cnd3d8fZ2fm+5O7h4YGXl5fR6+4AcXFx9dabgaqjk8nJyXz44YeN6pXaqVMnIiIiyMvLa9RfAc0tKCiIgQMHsn//fqusgClJlqbFJ3d9hcja1q2N7c6k5+/vT2BgYL31ZgBeffVVTp8+zf/+9z+j49bT91ytrKxk9erVJt/HEsTFxdGpUyfWrl3L8ePHmzscSbJqLT65Q+0VIqHqacqKigqTEo1Op6u33gxUtdLT6XS89957jap13qpVK3r37s2PP/5o0l8alkIIwdChQ/H09GTRokVcu3atuUOSJKslkzs/d2a697RGmzZtUKvVRq+7g2H1ZoQQvPbaaxw5coRFixYZPUZNvXr1wsvLixUrVlh1UwwHBwcyMjK4fv06ixcvttqNYklqbjK5U3UcEu7fVLW3tyckJMSk2bAh9WYAhg0bRseOHU1u5lFzvMGDB1NSUsLGjRtNvo8laN26NYMHD+bIkSN8++23zR2OJFklmdypmi16e3vXWgkyPDycc+fOmbRZaUi9GZVKxSuvvMIPP/zQ6GYcYWFhREZGsnHjRquvmR4TE0PXrl1Zv3690U8KS5Ikk3s1/abqvbPniIgIwLjuTHqG1JsBePzxxwkODja5mUdNAwcORKPRWHVhMahaskpLS8PHx4fFixdz5cqV5g5JkqyKTO536CtE3ruE4uvri4uLi8kblfp6Mz/99FOd19jb2/N///d/bNy4sdFLKq6urvTr148jR45YfUs7e3t7MjIyKC8vJysrS66/S5IRZHK/Q7+peu+RSFO7M+mFh4fj6elZ78YqwNNPP423t7fJzTxq0ul0+Pv7k5uba/UVF318fHjkkUcoLCwkLy+vucORJKshk/sdvr6+aDSaOtfdy8rKOHXqlNH3NbTejJOTE7/+9a9Zvnw5e/bsMXqcmlQqFWlpaVy9epV169Y16l6WoGvXrnTv3p2NGzdy8ODB5g5HkqyCTO53qNVq/P3960zuYHwpAr2YmBg0Gk2Ds/fnnnsOFxeXRjXz0AsMDCQ2NpZt27aZ9KVkaVJSUvDz8yM7O7ve00eSJFWRyb2GwMDA+ypEQtWs2pTuTHr6ejN79uypt9SAp6cnU6ZMYf78+WY5ITJgwACcnJxYvny5VW+uAtjZ2ZGRkcHt27dZtGhRvRvUkiTJ5H6XwMBAKioqau1PGh4eTlFRkcl1YHQ6XYP1ZgB+85vfoNFoeP/9900apyYHBwcGDhxIcXFxvccxrYWXlxdDhw6lqKjIquvYS9KDIJN7DXVtqsLP3ZnqKydQn4CAAIPqzQQEBDBx4kQ+//xzTp8+bdJYNUVFRREaGsratWtt4nH+Ll26EBcXx/fff1/vCSRJaulkcq/B3d0dJyenWtfdg4KCcHJyatRmpyH1ZqCqmcetW7eYNm2ayWPp6c+Ll5eXs2bNmkbfzxI8/PDDBAQEkJOTQ0lJSXOHI0kWSSb3GoQQBAUF1Zrc1Wq1QeUE6mNIvRmoenAqIyOD6dOnc+nSJZPGqsnb25uHHnqIH374gWPHjjX6fs1No9GQkZGBEIJFixZZdS0dSWoqMrnfIzAwkPPnz9e6tm5IOYH61Kw301A5g1dffZUrV64wffp0k8a6V2JiIh4eHixfvtwmNiM9PDwYPnw4p06dIjc3t7nDkSSLI5P7PfSdmWrr5+nh4UGHDh3YsWOHybPFHj16oChKgw2ho6OjSUlJYdq0aVy/ft2ksWqys7MjNTWV8+fPs3nz5kbfzxJ06NCBhx56iO3bt7N3797mDkeSLIpM7vfQJ/e6mk7rdDquX79Ofn6+Sff39PQ0qN4MVLXiO3fuHJ9//rlJY92rffv2dOzYkW+++cYsyz2WoH///gQHB/PVV18Z3e9WkmyZTO73cHBwwMvLq9Z1d6iqvNiqVasG183rY0i9GYA+ffrw0EMP8f777xvdDaouKSkpCCEaXYHSUqjVakaNGoWdnR2ZmZlm++8kSdZOJvda6DdVazuyKIRAp9NRVFRk8pOfhtab0TfzKCwsZP78+SaNdS93d3eSkpI4cOAABQUFZrlnc3Nzc2PkyJGcPXuWFStWNHc4kmQRZHKvRWBgINeuXavzVEx0dDR2dnZs3brVpPsbWm8GYPDgwURGRvL222+bbVaakJCAr68vq1atory83Cz3bG7h4eEkJiaye/dudu3a1dzhSFKzk8m9Fvp197qWZhwcHIiKimLfvn0mb3YaWm9GpVLx9ttvU1BQYJanVqFqKSMtLY3S0lI2bNhglntagqSkJEJDQ1mxYkWDX5qSZOtkcq+Fn58farW6zk1VgLi4OCoqKti9e7dJYxhabwZgyJAhjBo1ijfffJMDBw6YNN69QkJCiI6OZvPmzbWWW7BGKpWK9PR0HBwcyMzMtPpyx5LUGDK510JfIbK245B6fn5+hISEsH37dpObSBhabwbgo48+wsHBgcmTJ5utacXDDz+MVqu1+q5NNbm4uJCenk5JSQnLli2zmc8lScaSyb0OgYGBnDx5st7jinFxcVy8eNHkUsD6ejPbt29vMAn5+/vz/vvv88033zBr1iyTxruXk5MTycnJFBYWGvQFYy1CQ0Pp168f+/bts4mCaZJkCpnc6xAUFFRnhUi9jh074uLi0uhjkefPnzeoINlTTz1FUlISU6dONUtRMaha+w8KCmLNmjWUlZWZ5Z6WoHfv3kRERLBq1SqbqGcvScaSyb0ODW2qQtXyTY8ePTh06JDJBay6dOmCo6Njg0+sQtWa8owZMygrK+PXv/61SePdS19YrKyszKbK6AohGDFiBM7OzmRmZppcqlmSrJVM7nXw8PCos0JkTT169EClUpk8e9fXmykoKGiw3gxUPWX6+uuvk5mZydKlS00a816tW7cmPj6enTt31ruJbG2cnJwYNWoUpaWl5OTkyPV3qUWRyb0OQggCAwMbTO6urq506tSJ3bt3m3wO3diCZFOnTiUqKopnn33WoC8EQ/Tt2xdXV1eWLVtmtg1bSxAcHExycjIFBQVs2bKlucORpAdGJvd6BAYGcu7cuQaP1Ol0Om7cuGFy8Sp9vZkdO3YYVLHR3t6ezz77jJMnT/K73/3OpDHvpdVqSUlJ4cyZMyY/nGWpEhIS6NChA2vWrLGpv0wkqT4yuddD35mpodl7SEgIvr6+DXZZqk9sbKxB9Wb04uPjef7555k+fTrfffedSWPeq1OnTkRERJCXl2e2vwgsgRCCYcOG4ebmRmZmplmqbEqSpTMouQshUoQQ+4UQh4QQr9ZzXboQQhFCxJovxOYTEBAANJzc9fVmTp8+bfLMMCIiwqB6MzX95S9/ISgoiEmTJpnlgR0hBKmpqVRWVrJ69epG38+SODo6kpGRwbVr11iyZIlcf5dsXoPJXQihBv4NpAKdgceEEJ1ruc4VeAGwmYVNR0fHeitE1tS1a1e0Wu0DqTej5+rqyscff0x+fj7vvvuuSePeq1WrVvTu3Zsff/yRw4cPm+WeliIgIIBBgwZx8OBBNm3a1NzhSFKTMmTmHgccUhTliKIo5cB8YFgt1/0/4F3Aps6c6TdVG5rp2dvbEx0dTX5+PlevXjVpLEPrzdSUlpbGmDFjeOutt8zWMLpXr160atWKFStW2FwLu9jYWLp06cK6detsouWgJNXFkOQeCJyo8XPRnd9VE0J0B4IVRVle342EEJOFENuFENvPnTtndLDNITAwkKtXrxq0Bq3T6aisrGTnzp0mjWVMvZmapk2bhrOzM5MmTTLLSReNRkNaWholJSVs3Lix0fezJEIIhgwZQqtWrcjKyjL5i1iSLF2jN1SFECrgQ+Dlhq5VFGWGoiixiqLE+vj4NHboB6Khzkw1eXl5ER4e/sDqzej5+fnx4YcfsmnTJmbMmGHSuPcKCwsjMjKSjRs32lyHI61WS0ZGBjdu3GDx4sU2dfRTkvQMSe7FQHCNn4Pu/E7PFYgE1gshjgEJwFJb2VRt3bo1arXaoHV3qErOV65cMbkRhjH1ZmqaOHEiAwYM4JVXXjE41oYMHDgQe3t7MjMzbabuu56fnx+DBw/m6NGjNlX2WJL0DEnu24B2Qoi2Qgh7YAxQ/WikoiiliqJ4K4oSqihKKPA9MFRRlIafp7cC+gqRhibMdu3a4e7ubpZ6M8asCQsh+OSTTygvL+e5554zeeyaXF1dSU9P59y5c3z11Vc2d8IkOjqabt268c0339jc5rEkNZjcFUWpAJ4DcoGfgIWKovwohHhTCDG0qQO0BIGBgZw6dcqgP99VKhWxsbEcO3bM5Drp+nozxn5BRERE8Kc//YklS5awePFik8a+V3h4eHWFRVt7wlMIweDBg/Hx8WHx4sU2dbZfkgxac1cUZYWiKO0VRQlXFOWtO797Q1GU+4qbKIrS11Zm7XqBgYHcunXL4GTdvXt31Gr1A6s3U9NLL71EdHQ0zz33HJcuXTJp/Hv17t2bjh07snr1agoLC81yT0thb29PRkYGt27dIisrS66/SzZDPqFqAGM2VaGqYJX+1IupDxcZW29Gz87OjpkzZ3LmzBleeeUVk8a+lxCC4cOH06pVKzIzM7ly5YpZ7mspfHx8GDJkCMePH2fdunXNHY4kmYVM7gbw9PQ0qEJkTTqdjvLycpObYHh6etKuXTuD683U1KNHD1588UVmzJhhts1CrVbL6NGjKS8vZ+HChUbHZOmioqLo0aMHmzZtMlsrQ0lqTjK5G8DQCpE1BQYGEhAQ0Kh6Mzqdzqh6MzW9+eabhIaGMnnyZLPVMvfx8WHYsGEUFRWRm5trlntakpSUFFq3bk12djYXLlxo7nAkqVFkcjeQoRUiazKmy1JtTKk3o+fs7Mynn37K/v37efvtt00avzZdunShZ8+ebNu2zaZa80HVXkdGRgYqlYo5c+bIDVbJqsnkbiD9unt9TbPvFRkZadKpFz1T6s3UNHDgQMaNG8c777zDvn37TIqhNsnJyYSGhrJs2TKztfuzFK1atWLs2LGUlZUxZ84cWUFSsloyuRtIn9yNOQ+t0Wjo3r07+/fvp7S01KRxo6Ojja43U9OHH36Iu7s7kyZNMts6uUqlYtSoUTg6OrJgwQKb6r0KVQ+SPfbYY1y8eJG5c+eapeKmJD1oMrkbyNHRkc6dO7Nt2zajZnOxsVUP6hrSI7U2+pM3u3fvNmn27uPjw7Rp0/j+++/5+OOPTYqhNs7Ozjz66KNcvnyZxYsX29wDTqGhoWRkZHDq1CkWLFhgcwXUJNsnk7sRkpKSKC8vZ/PmzQa/x8PDg/bt27Nz506TE0RycjIODg5kZWWZ1Mpv7NixDBo0iNdee40TJ040/AYDBQUFkZqayqFDh1i/fr3Z7mspOnTowLBhwzh69Kg8Ay9ZHZncjeDr60tkZCRbtmzh2rVrBr9Pp9Nx/fp18vPzTRrX2dmZESNGcO7cOZOaaAgh+Pjjj6msrOTZZ5816yy7R48eREdHs2HDBps8QtitWzdSUlIoKCiwyRIMku2Syd1ISUlJVFRUGNXaLiwsDC8vr0bVmwkPD6dnz55s377dpKJkbdu25f/9v//HsmXLyMzMNDmOe+kf4ff392fx4sWUlJSY7d6WIj4+nqSkJHbv3k1ubq5M8JJVkMndSN7e3kRFRbF161aDa4HrT70UFRUZddrmXgMGDMDf35+lS5eadEzv17/+NT169OD55583axK2s7Pj0UcfRaVSsWDBApurIAlVX+pxcXFs2bKFb7/9trnDkaQGyeRugsTERG7fvm1Uq7bo6Gjs7OwaNXtXq9Wkp6dTUVFBdna20WvAGo2GmTNncuHCBaZOnWpyHLXx8PAgPT2ds2fPsmzZMpub3QohSElJoVu3buTl5ZncTlGSHhSZ3E3g5eVFt27d2L59u8F1VhwcHOjatSv79u1r1NlpLy8vUlNTOXbsmEl9QKOjo3n55ZeZNWsWeXl5JsdRG30Fyb1799pk8hNCMHToUDp06MDKlSvZs2dPc4ckSXWSyd1E+tm7MW3odDodFRUV7N69u1FjR0dH06VLF/Ly8gwuZlbTH//4R8LDw5k8ebLZz6j36dOHDh06sHr1ao4fP27We1sC/Rn/0NBQlixZYpObyJJtkMndRJ6enkRHR7Njxw6D17/9/PwICQlh27ZtjTpWJ4TgkUcewc3NjaysLKMfsnFycuLTTz/l0KFDvPnmmybHUVdsw4cPx8PDwyYrSELV8taYMWPw9/cnMzNTNtqWLJJM7o2QmJiIoihGbbDFxcVx6dIlDh061KixHRwcGDlyJKWlpSxfXm9f8loNGDCAJ554gvfff9/sNWIcHBwYPXo0N2/eZNGiRTZXQRKqqmSOHTsWDw8PvvzyS06dOtXcIUnSXWRybwQPDw9iYmLYuXOnweUFOnbsiIuLS6M2VvVCQkJISkpi7969JiXov/3tb7Rq1Yqnn37a7AnY19eXoUOHcvz4cZPO5lsDJycnxo8fj6OjI1988YXNNRKXrJtM7o3Up08fhBAG101Xq9X06NGDQ4cOmeU4Yp8+fQgJCWHFihVG38/Ly4uPPvqI7du389FHHzU6lntFRkaSkJDA1q1b2bt3r9nvbwnc3NwYP348QgjmzJljcg0hSTI3mdwbyd3dne7du7N7926D29r16NEDlUplltm7SqVi5MiRqFQqFi9ebPQMfPTo0QwePJg//OEPTbJ2nJycTJs2bVi6dKlJtXGsgZeXF+PGjePmzZvMmTPHqKeXJampyORuBr179zZq9u7q6kqnTp3YvXu3SbVi7uXu7s6QIUMoLi42+nijvjSBEIJnnnnG7OfT1Wr1XRUkzdU4xNK0bt2axx9/nNLSUr744gub/ZyS9ZDJ3Qzc3NyIjY1l9+7dBi+N6HQ6bty4Ybblis6dOxMTE8OmTZuMbg4SEhLC22+/zapVq/jyyy/NEk9NLi4uZGRkUFpaSnZ2ts094KQXEhLCo48+ytmzZ5k/f75ZvrglyVQyuZtJ7969UavVBs/eQ0JC8PX1bVQbvnulpKTg7e3N4sWLjX5Q6le/+hXx8fG88MILTbIxGBwcTEpKCgcOHDBbX1dL1K5dO0aMGEFhYaHNnhSSrINM7mbi4uJCbGwse/bsMaj/phACnU7H6dOnzVaG197envT0dMrKyli6dKlRXxpqtZrPPvuMS5cu8fLLL5slnnvFxsbStWtX1q9fz8GDB5tkDEsQGRlJWloaBw4cICcnx2b/UpEsm0zuZtSrVy80Gg3ffPONQdd37doVrVZrlo1VvdatW5OcnMz+/fuNbhASFRXFb3/7W2bPns2aNWvMFpOe/uErPz8/Fi9ezMWLF80+hqWIjY2lf//+7N27l5UrV8oELz1wMrmbkYuLCzqdjr1793Lu3LkGr7e3tyc6Opr8/HyDK0waIj4+noiICFavXs3Zs2eNeu/rr79O+/bt+eUvf9kkpz7s7OwYPXo0AAsWLLDpdenevXtXNxO3xWYmkmWTyd3MevXqhb29vcHryjqdjsrKSnbs2GG2GIQQDBs2DK1Wy6JFi4xKoA4ODsyYMYOjR4/ypz/9yWwx1eTp6cnIkSM5c+aMTVaQ1BNC8PDDDxMTE8OGDRv4/vvvmzskqQWRyd3MnJyciIuLY9++fQbNmr28vAgPD2fHjh1m3XxzcXFh+PDhJnVvSkpK4umnn+bDDz9k586dZouppnbt2tG3b1/27Nlj1mUpS6NfiurUqRO5ubmNLhonSYaSyb0JPPTQQ9jb2xu89q7T6bhy5Qr79+83axwREREkJCSwfft2o+/93nvv4evry9NPP91kzaETExNp3749ubm5Zu3tamn0D5qFhYWxdOlSkzppSZKxZHJvAo6OjiQkJJCfn8/p06cbvL5du3a4u7s3yQx2wIABtG7dmpycHKO6N3l6evLPf/6TXbt2MWnSpCY50ieEYMSIEbi7u7Nw4UKz7jtYGo1Gw+jRowkMDGTRokUcOXKkuUOSbJxM7k0kISEBrVZr0OxdpVIRGxvLsWPHjN4AbYhGo6nu3rRkyRKjSg2PGjWKP/7xj/z3v/9l7NixTbL5qa8geePGDTIzM236XLi9vT2PP/44Xl5ezJ8/n+Li4uYOSbJhMrk3EUdHR3r27ElBQYFB5WC7d++OWq1uktm7t7c3KSkpHD161KjG3gB/+tOfeO+991iwYAGjRo1qksfq/fz8qitINsURTEvi6OjIuHHjcHFxYe7cuQadqpIkU8jk3oTi4+NxcHAw6Bick5MTkZGR7Nmzx+jmG4aIiYmhc+fO5OXlGT1jnDp1Kv/6179YunQpQ4cObVSbwLpERUVVN6C21QqSeq6urowfPx61Ws2cOXMMLjgnScaQyb0JOTg48NBDD3HgwAGDEqpOp6O8vNzszTOgan17yJAhuLq6mtS96Ve/+hWzZs1i7dq1pKamNkmHpYEDBxISEsJXX31l9uUpS+Pp6cn48eO5desWc+bMsen9Bql5yOTexOLi4nB0dDRo9h4YGEhAQIBZ683UpO/edOnSJVasWGH0+5988knmzp3Lpk2bSE5ONvsTpvoKklqt1qYrSOr5+voyduxYrly5IitJSmZnUHIXQqQIIfYLIQ4JIV6t5fWXhBD5Qog9Qoi1Qog25g/VOmm1Wh566CEOHTpk0HG/uLg4zp8/b3RlR0OFhISQmJjInj172LNnj9HvHzNmDFlZWezevZt+/fqZfYbt6upKRkYGly5dYuHChTaf8IKCghg9ejTnzp1j3rx5lJeXN3dIko1oMLkLIdTAv4FUoDPwmBCi8z2X7QJiFUXpCiwC3jN3oNYsLi4OJycng2bvXbp0wcnJqUkf7ElMTCQ4OJjly5ebNPseNmwYX331FQcOHCApKYmTJ0+aNb6QkBCGDh1KYWEhn332mc1vOoaHh5Oenk5RUZFcg5fMxpCZexxwSFGUI4qilAPzgWE1L1AUJU9RFP0u2/dAkHnDtG729vb06tWLI0eOUFhYWO+1Go2GmJgY9u/f32Qt2/QP1QghyMrKMun44cCBA1m1ahVFRUUkJiY2+LmM1a1bNyZOnMjNmzeZOXOmzT/407lzZ0aNGsW5c+f45JNP2LdvX3OHJFk5Q5J7IFBzPaHozu/q8hSwsrYXhBCThRDbhRDbbX02di+dToezs7NBs/fY2FgAo6s6GsPDw6O6e5OpRa0SExP5+uuvuXDhAn369DF7Gd+QkBAmT56Mt7c3CxYsYP369TZbhwaqEvwvf/lLfHx8yMrKYunSpXKZRjKZWTdUhRDjgFjg/dpeVxRlhqIosYqixPr4+JhzaItnZ2dH7969OXbsWIO9Sj08PGjfvj07d+5sskf/oWoJKCYmho0bN5q8xh8fH09eXh5lZWUkJiby448/mjVGNzc3nnzySbp168Y333zDggULmuSoqKXw9PTkiSeeoE+fPuzatYsZM2YY9JSzJN3LkOReDATX+Dnozu/uIoRIBn4PDFUUxXb/39cIPXr0wNXV1aAZqE6n4/r16+Tn5zdpTCkpKXh5eZGdnW3y+fXo6Gg2bNiAEIKkpCSzFxvTaDQMGzasupPTzJkzDWqIYq3UajX9+/dnwoQJlJeXM3PmTL7//nub/qtFMj9Dkvs2oJ0Qoq0Qwh4YAyyteYEQIgb4lKrEbtsHlBtBP3svLCxscKYcFhaGl5dXk1dM1HdvunbtmtHdm2rq1KkT3377LS4uLvTv35/NmzebNU4hBPHx8YwfP55r167x2Wef2XQ3J4C2bdsyZcoUwsPDyc3N5csvv2ySGvuSbWowuSuKUgE8B+QCPwELFUX5UQjxphBi6J3L3gdcgEwhxG4hxNI6btfide/eHTc3twZn70IIYmNjKSoqMvtplHv5+/ub3L2ppvDwcDZs2ICPjw8PP/wweXl5ZoyyStu2bZk8eTKenp7MmzePjRs32vSM1snJiTFjxpCamsqRI0f45JNPZNExySAGrbkrirJCUZT2iqKEK4ry1p3fvaEoytI7/05WFMVPUZToO/8bWv8dWy6NRkOfPn04ceIEhw8frvfa6Oho7OzsHki984SEBMLDw03q3lRTSEgIGzZsoE2bNgwePJhVq1aZMcoqHh4e/OIXvyAyMpK1a9eSlZVl0xuPQgji4uKYNGkSDg4OzJkzh6+//tqmi6xJjSefUG0GMTExuLu7Nzh7d3BwoGvXruzbt69J6rnUJIRg+PDhaLVasrKyGlUB0t/fn2+++YZOnToxdOhQsrOzzRhpFTs7O0aOHElycjL5+fnMmjXLpnuyQlWBtcmTJ9O9e3c2bdrE559/bvOfWTKdTO7NQK1Wk5iYSHFxcYPrxjqdjoqKCnbt2tXkcem7N509e7bR1Rm9vb1Zt24dsbGxZGRkMG/ePDNF+TMhBL169eLxxx+ntLSUzz77zOaXLOzs7BgyZAgZGRlcuHCBTz75xOYLrUmmkcm9mXTr1g0PD48GZ+9+fn60adOGzZs3c+bMmSaPS9+9adu2bY3uDOXh4cHq1avp06cP48aNY+bMmWaK8m4RERFMmjQJFxcXvvjiCzZv3mzT6/Dw85l4Pz8/Fi9ezJIlS2z6iKhkPIw2lpQAABjRSURBVJncm4larSYpKYlTp05x4MCBeq9NTU1FpVLxn//854E8qVmze1Njqz+6uLiwYsUKUlJSmDRpEh999JGZorxbq1ateOqpp+jYsSOrV69myZIlTdJcxJJ4eHjwxBNPVNcKmjFjRpNvvkvWQyb3ZtS1a1datWpFXl5eg7P3SZMm4ePjw4IFC/j222+bdGZas3tTdnZ2o8dydHQkOzubESNG8MILL/DOO++YKdK7abVaMjIy6NevH3v27OHzzz9vshIOlkKlUtGvXz8mTpxIRUUF//nPf/juu+9s/i8XqWEyuTcjlUpFYmIiZ86caXBG7urqyhNPPEFUVBTr1q0jOzu7SWemNbs3zZs3z6j+q7XRarUsXLiQsWPH8rvf/Y7XX3+9SRKQEILExETGjBlDSUkJM2bMMHvdG0vUpk0bpkyZQvv27VmzZg1z586VNeJbONFc3/CxsbFKU9ZOsRaVlZVMnz4dtVrNlClTEELUe72iKGzcuJF169YREBDAmDFjcHV1bZLYFEVh27ZtrFmzBrVaTUpKCt26dWswxvrcvn2bKVOmMHPmTH7zm9/wwQcfNOp+9Tl//jzz58/n4sWLDBo0CJ1O12RjWQpFUdixYwe5ublotVqGDx9OREREc4clmZEQYoeiKLENXieTe/Pbu3cvixcvZtSoUXTp0sWg9xQUFLB48eLqBtOBgfXVcmuckpIScnJyOH78OO3atavu6GQqRVF48cUX+eijj/jlL3/J9OnTUama5o/IGzdukJ2dzYEDB4iOjiYtLQ2NRtMkY1mSs2fPkpWVxdmzZ+nZsycDBgxArVY3d1iSGcjkbkUqKyv55JNPUBSFZ555xuBEd+bMmepH0ocNG0ZkZGSTxagoClu2bGHt2rVoNBpSUlLo2rWryTNhRVH4/e9/zzvvvMP48eOZNWtWkyVdRVHIy8vj22+/JTAwkNGjRzfZXzuW5NatW6xevZrt27fj7+9Peno6Xl5ezR2W1EgyuVuZH3/8kUWLFjFy5EiioqIMft+1a9dYuHAhx48fp0+fPvTr169Jlx4uXLhATk4OJ06coEOHDjzyyCO4uLiYfL+33nqLP/zhD4waNYq5c+dib29vxmjvlp+fz5IlS9BqtTz66KMEBwc3/CYbUFBQQE5ODrdv3yYtLa1RX8pS85PJ3cooisInn3zC7du3efbZZ41aprh9+zbLly9n165ddOzYkREjRjRpkqysrGTLli2sW7cOOzs7UlNTiYyMNDlhTJs2jd/85jekpaWxaNEiHBwczBzxz86ePcv8+fMpLS0lLS2N7t27N9lYlqS0tJTs7GwKCwuJiooiLS0NrVbb3GFJJpDJ3Qr99NNPLFy4kOHDh9OtWzej3qtfNlm9ejW+vr6MGTMGDw+PJoq0yvnz58nJyaGoqIiOHTuSlpZm8ix+xowZTJkyhf79+5OTk4Ozs7OZo/1ZWVkZWVlZHD58mNjYWFJSUlrEenRlZSUbN25k/fr1eHh4kJ6e3qR7NVLTkMndCimKwowZM7h58ybPPfecSZuMhw4dYtGiRajVakaPHk1ISEgTRPqzyspKNm/eTF5eHvb29gwePJguXbqYNIv/4osvmDhxIj179mT58uW4u7s3QcRVKisrWbt2Ld999x0hISFkZGQ0annJmhw/fpzFixdz5coV+vXrR69eveQyjRWRyd1K7d+/n/nz5zN06FBiYmJMusf58+f58ssvuXTpEo888ojJ9zHGuXPnyMnJobi4mE6dOpGWlmbS7HvRokU89thjREdHs2rVqibfANy3bx85OTk4OTkxevRoAgICmnQ8S1FWVsayZcvIz88nLCyM4cOHt4hNZlsgk7uVUhSFzz77jLKyMp577jmTlwvKyspYtGgRR44cISEhgYcffrjJjhvqVVZW8t1337F+/Xq0Wm31LN5Yy5cvJz09HW9vb6ZOncrTTz/dpMs0p0+fZv78+Vy9epUhQ4YYvSRmrRRFYdeuXaxcuRKNRkOPHj3o0aMHnp6ezR2aVA+Z3K3YwYMHmTdvHkOGDGnUhl9lZSW5ubls3bqV8PBwRo0a1aSblXpnz55lyZIlnDp1ii5dujB48GCcnJyMusf333/Pq6/+//bOPbat677jnyOKT1GyZMt6xLJkOXH8iCVbDzpFOwzt2szRbKcZujYxumEoAhSrlzTu0K5J/um2P+oURdMGDVC0aVokXjA3SbMo8IJsKdzUVrolpB62/K5sPSyJFEXrwafE19kfEi9EWbYlmRQt8nyAA957dXXu757L+z2/+7u/c/gMf/jDH1i3bh1PP/00Tz75ZNqEJxgM8uabb9LX18eDDz7I5z73uZx54Tg6OsqJEye4dOkSUkq2bNmCzWbj3nvvTbtDoFg6StxXMVJKXnnlFfx+P0899dQdv+xrb2/nvffeo6SkhIMHD65IrnM8Huejjz7iww8/xGw2s2/fPrZv377kev74xz9y5MgRjh8/jtVq5Rvf+Abf+ta3qKysTLnNsViMDz74gI8//hiDwUB9fT179uwhV37M3ev10t7eTkdHB36/n+LiYpqbm2loaFhy56xIH0rcVzk9PT28/vrr7Nu3j+bm217H29LX18cbb7yBlJIvf/nLbN68OQVW3p6RkRFaW1txOp3s3LmTlpaWZQlFd3c3zz//PMeOHUOv1/O1r32N73znO2k5j6GhIex2O2fPniUWi7Fp0yZsNhvbtm3LCU82Fotx8eJFHA4HfX196HQ6HnjgAWw2Gxs2bFAvXzOMEvdVjpRSm9XwqaeeSsnozfHxcY4dO8bo6CgPP/zwis21EovFaGtr4+TJk5jNZvbv38+2bduWVdeVK1f44Q9/yK9//WtisRiPP/44zzzzTFpG5wYCATo7O3E4HExOTlJUVERTUxONjY05k1njdrtxOBycPn2acDhMRUUFNpuNnTt3pnUsheLmKHHPAq5evcrRo0dpaWlhz549Kalzenqat99+m8uXL9PU1ERLS8uK5Xi7XC5aW1txuVzU1dXR0tKC2WxeVl3Dw8P8+Mc/5mc/+xmBQIADBw7w3HPP8alPfSrFVs+EmC5fvozdbufq1avk5eVpnmxVVVVOeLLT09N0d3djt9txu90YjUZ2795Nc3MzpaWlmTYvp1DingVIKXn11Ve5fv063/zmN9Hr9SmpNx6Pc+LECT766CNqamr4yle+smIx1VgsxqlTpzh16hQWi4X9+/ezdevWZdc3NjbGSy+9xIsvvsjY2Bif/exnefbZZ3nooYfSIroejwe73c7p06eZnp7WPNm6urqUXZ+7GSkl165dw263c/78eeLxOLW1tdhsNrZu3ZoTYatMo8Q9S+jr6+PVV19l7969KfdKz5w5w7vvvkthYSEHDx6krKwspfXfCqfTyTvvvIPb7WbXrl3s3bt32V48gN/v5+WXX+ZHP/oRQ0NDNDU18dxzz/Hoo4+mRXDC4TBnzpzRPFmTyURDQwPNzc2sXbs25ce7G/H7/VrYyuv1UlhYqIWtVM58+lDinkW89tpr9Pf3s2PHDpqbm6murk6ZVzo4OMhvfvMbwuEwX/rSl7j//vtTUu9iiMVinDx5klOnTmG1Wtm/f/8dH396epqjR4/ygx/8gJ6eHrZt28Z3v/tdvvrVr6bFs5ZSMjAwwCeffMKFCxeSUgnvu+++nAjZxONx/vSnP2G327ly5Qp5eXls27YNm81GTU1NTrTBSqLEPYvw+Xy0tbVpoYCysjKamprYtWtXSnKxvV4vx44dw+l08oUvfIFPf/rTK3pDDg8P09raitvtZvfu3ezdu/eO8/FjsRhvvfUWR44c4fTp01RXV/Ptb3+bJ554Im0hKJ/PR3t7O+3t7fj9fkpKSrRUwjt5KllNjI2N4XA46OzsZGpqitLSUmw2G/X19SsyxiIXUOKehYTDYc6ePYvD4cDpdKLX66mrq6O5ufmO874jkQitra2cO3eO+vp6Dhw4sKI/ahGNRjl58iRtbW1YrVYeeeSRlPyCkJSS999/n+9///u0tbWxfv16Dh8+zKFDh9I2sVosFuPChQvY7XYGBgbIz8+nrq4Om82Wlvz8u5FIJMK5c+ew2+0MDw+j1+upr6/HZrNRXl6eafNWNUrcs5yhoSEcDgdnz54lGo2yYcMGmpubeeCBB5YdfpBScvLkST788EOqqqp47LHHVjzlb2hoiNbWVkZHR7n//vvZsmULmzdvpqSk5I6fJtra2jhy5AjvvfceRUVFHDp0iMOHD6dVbFwuF3a7ne7ubiKRCFVVVezZs4cdO3bkxEyUcON3dePGjdhsNrZv354Tv4qVarJX3A8fhq6u1Bu0SonF4wT8fnw+H5FIhLy8PKxWK4WFhcsW+UAwiMfjIS8vj7KyMowrnM8spWRichK/308sGgUgPz8fk9mM2WTCZDaju4OXpH6/n4GBAdyjo+QJQUVlJdUbN6Y1bBCPx/H7/Xh9PqKRCDqdDmthIYWFheTniMjHZtvAN9sGeTodhVYrFosFg8GQW7H53bvhJz9Z1r8uVtxVt7nK0eXlUVRURGFREVNTU/h8Prw+H16vF5PZTGFhIRazeUk3ToHFgr6ighG3G5fTicViwWQyYTKZyNfrSfctKISgpLiY4uJiopEIoakppkIhgoEAfp8PAIPBoIm90WQibwnnZ7Va2bFjB5tCIa4NDOB0OnEOD1NWXk51dTUFaYjJ5829TqEQXp+PyYkJJicnsVgsFBYWYjKZ0t62mUSXl8eaoiKKZtvA5/Mx6fUyOTkJQmDQ6zEYjRgNBgxGIwa9PrcEP8WsPs9dcVt8Ph+dnZ20t7fj9XqxWq00NjbS1NREUVHRouvx+/387ne/48qVK/j9fmBGGGtqaqipqWHTpk2Ulpau2A0Yj8cZHh7m6tWrXL16lWvXrhGPx9HpdFRXV1NbW8vmzZuprKxcUvrj4OAgL7zwAj//+c8JBoMcOHCAz3/+8zQ2NrJr164ltdlSGB8f114+hkIh1q9fz/bt26moqKC8vDwloai7HZ/Px+DgIMPDw1qZmpoCQKfTUV5ezj333KOV9evX53wuffaGZRSLJpGi5nA46OnpQQjB1q1baWpq4t577120cEgpGRsbo6+vj/7+fvr6+vDNetAWi0UT+pqaGsrKylZMkMLhMAMDA1y5coXe3l5GRkYAMJlMmtAvJV7v8Xj46U9/yssvv4zT6dS2b9myhYaGBhobG2loaKChoSGlk4klXj46HA6Gh4dJ3JMGg4Hy8nKtVFRUUFZWltXD/qWUTExMMDw8zNDQEE6nk+HhYcLhMDATnqusrKSyslIT/HXr1uWU4CtxVyQxPj5Oe3s7nZ2dBINBSkpKaGpqWtaMf1JKxsfH6e/v18R+cnISALPZnOTZl5eXr5jY+/1+ent7Nc/e6/UCUFxcrIl9bW3touaGdzqddHZ20tHRQWdnJ52dnfT29mp/r6qqShL8xsbGlExFEIlEcLvdjIyM4HK5GBkZYWRkhOnpaW2ftWvXat59QvSLioqy1suXUnL9+nXNs3c6nTidTiKRCDDTCc4X/LVr12ZteyhxVyxINBrlwoULOBwOBgYG0Ol02uCojRs3LvuGmJiYSPLsJyYmgBkver7Yr4SXlXjaSAh9b2+vJpAVFRWaV19dXb3oF8/j4+N0dXVpgt/R0cGlS5eIx+MArFu3Lsm7b2xs5L777rvj85VSMjk5mST2LpeL8fFxbR+TyZQk9uXl5ZSVlWVtNko8Hsfj8SSFc1wuF7FYDACj0ZgUzrnnnntYs2ZNVgi+EnfFbZk/419ZWRnNzc3U19ff8eCoyclJTez7+/sZGxsDZm666upqTeyXGh9fLnPj9b29vQwMDKQkXh8IBDhz5ozm3Xd0dHD27FktjGC1Wtm9e3eSl79jx46UjJadnp7G7XYnif7IyIjm0QohKC0tvUH0rVZrVojcfGKxGKOjo0mCPzIyonW+ZrM5SewTbbHa5gRS4q5YNOFwmO7ubhwOBy6XC4PBoA2OqqioSMkxvF6v5tX39/dz/fp1YOaRer7Yr0T+dyJen/Ds58bra2pqWLNmDQUFBVitVgoKCpKWbycG4XCY8+fPJ4V1urq6CAQC2jnX1dUlCX59fX1KRs4mnljmhnVcLpcWooKZ9yTzwzqlpaVZmXcfjUZxu90MDQ1pIR23281c3dPr9UnXeX6Z+zfzEjPP0oESd8WSkVJqA07OnTunDY4qKyvDYrFgsVgwm803LJtMpiV73z6fLylm7/F4gJkbLSH25eXlWgpmoujTlB4XCAS0eP21a9fw+XxJce65GAyGG276m3UERqMRIQSxWIyenp4kwe/o6NCeaBLzsdTW1lJSUnLLUlxcTElJCRaLZdFtEQqFbojju91uLYwBMx2b2WzWPucvzy+Jv602zzcSieByufB4PAQCgaTi9/sJBAIEg0EW0kYhxA3X3GKx3LRzSEdYLKXiLoR4GHgR0AG/lFI+P+/vRuA1oAm4Djwmpey7VZ1K3O9uQqEQXV1ddHd34/P5CAaD2uPtQiSEfr74L9QZJJbneoqBQCDJs3e73QseRwhxg+CbTCaMRuOC2+eXpQyWiUajSTf8/OX5YrAQOp1uwY4gUYLBIL29vVy8eJHOzk4GBwcZHx9nYjYH/lb3p16vv20HcLNitVqRUuLxeBgZGcHj8RAKhZiamiIUCt2wfCs78vPzF+wE5m+bv240Gu/aLJd4PE4oFFrwWi/UIURnB9vNx2g0Lnjtt27duuypKFIm7kIIHXAZeAgYBOzAQSnl+Tn7HALqpZT/IIR4HPhrKeVjt6pXifvqQkpJOBwmGAwSDAYJhUI3LM//DAaDN/3Sw8wX/2adQX5+PpFIhGg0qpXEeiQS0Uo4HNY+w+HwLY8HM52DwWDQOoO5nUJCgBJPInl5eQgh0Ol0CCGSts1fhplwzPT0NNPT00xNTWllrlAm2mehjjJhW15eHvn5+Vr9ifaPx+PEYrGk9giHw0nHCgaDWmcTiUSIxWJJJRqNEovFkFJqbV1QUIDFYkGv16PT6dDpdOTn599QDAZD0vr8NkhoiZQSKaV2nFuh1+u1unQ6nVZnYjlhz9zlhdbz8/OT1udum/s5f1ui3eeXpWyHmaeBREc4975IXIu51yUUCnHgwAEaGxtv2Ta3+A6nbITqHqBHSnl1tuJjwBeB83P2+SLwL7PLbwEvCSGEzFTMR5FyhBAYjUaMRiMlJSWL/r9IJHLTzmD+ssfjIRgMai8j04GUUhPguXHou4GEbYtFp9NpHWImiMfjt3yaWwyJTnohliofmY6FL4Xjx48vW9wXy2LEfQNwbc76IPDgzfaRUkaFEJPAOsAzdychxNeBr8+u+oUQl5ZjNFA6v+4cQJ1zbqDOOTco/d73vrfcc65ZzE4rmgQrpfwF8Is7rUcI4VjMY0k2oc45N1DnnBusxDkv5m3GELBxznrV7LYF9xFC5ANrmHmxqlAoFIoMsBhxtwNbhBC1QggD8Djw7rx93gX+fnb5b4ATKt6uUCgUmeO2YZnZGPqTwH8zkwr5KynlOSHEvwEOKeW7wCvAUSFEDzDGTAeQTu44tLMKUeecG6hzzg3Sfs4ZG8SkUCgUivRxd44gUCgUCsUdocRdoVAospBVJ+5CiIeFEJeEED1CiGcybU+6EUL8SgjhFkKczbQtK4UQYqMQ4vdCiPNCiHNCiKczbVO6EUKYhBCfCCFOz57zv2bappVACKETQnQKIY5n2paVQAjRJ4ToFkJ0CSHSOkR/VcXcFzMVQrYhhPhzwA+8JqXcmWl7VgIhRCVQKaXsEEIUAu3Ao1l+nQVQIKX0CyH0QBvwtJTy/zJsWloRQvwT0AwUSSn3Z9qedCOE6AOapZRpH7S12jx3bSoEKWUYSEyFkLVIKU8yk4GUM0gpnVLKjtllH3CBmVHQWYucwT+7qp8tq8fzWgZCiCpgH/DLTNuSjaw2cV9oKoSsvulzHSHEJqAB+DizlqSf2RBFF+AGPpBSZvs5/wT4Z+DOJqhZXUjgf4QQ7bPTsaSN1SbuihxCCGEFfgscllLeXbN8pQEpZUxKuZuZUeB7hBBZG4YTQuwH3FLK9kzbssL8mZSyEWgB/nE27JoWVpu4L2YqBEUWMBt3/i3wupTy7Uzbs5JIKSeA3wMPZ9qWNPIZ4JHZGPQx4C+EEP+eWZPSj5RyaPbTDfwnM6HmtLDaxH0xUyEoVjmzLxdfAS5IKV/ItD0rgRBivRCieHbZzEzSwMXMWpU+pJTPSimrpJSbmLmPT0gp/zbDZqUVIUTBbIIAQogC4C+BtGXBrSpxl1JGgcRUCBeAN6SU5zJrVXoRQvwH8L/AViHEoBDiiUzbtAJ8Bvg7Zry5rtnyV5k2Ks1UAr8XQpxhxon5QEqZE+mBOUQ50CaEOA18AvyXlPL9dB1sVaVCKhQKhWJxrCrPXaFQKBSLQ4m7QqFQZCFK3BUKhSILUeKuUCgUWYgSd4VCochClLgrFApFFqLEXaFQKLKQ/wdGr1gR2H62XQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c8b5b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "testmus = np.linspace(0, 5, 11)\n",
    "results = []\n",
    "for mu in testmus:\n",
    "    results.append(runOnePoint(mu, data, pdf, init_pars, par_bounds))\n",
    "\n",
    "obs = [1.0 / r[-2:][0] for r in results]\n",
    "exp = [[1.0 / r[-2:][1][i] for r in results] for i in range(5)]\n",
    "\n",
    "\n",
    "def plot_results(testmus, cls_obs, cls_exp, test_size=0.05):\n",
    "    plt.plot(testmus, cls_obs, c='k')\n",
    "    for i, c in zip(range(5), ['grey', 'grey', 'grey', 'grey', 'grey']):\n",
    "        plt.plot(testmus, cls_exp[i], c=c)\n",
    "    plt.plot(testmus, [test_size] * len(testmus), c='r')\n",
    "    plt.ylim(0, 1)\n",
    "\n",
    "\n",
    "plot_results(testmus, obs, exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}