{ "cells": [ { "cell_type": "markdown", "id": "f1ce3bc5", "metadata": {}, "source": [ "# Maximum-likelihood estimation of $\\nu$ and $p$\n", "\n", "This example fits the two parameters of the Housworth–Stahl model to a set of crossover positions:\n", "\n", "- $\\nu$, the interference parameter ($\\nu = 1$ is no interference);\n", "- $p$, the fraction of crossovers escaping interference.\n", "\n", "It uses simulated data so the estimates can be compared with the truth. To use your own data,\n", "replace the simulation step with a list holding one array of crossover positions (in Morgans) per\n", "meiosis, and the chromosome length(s) `L`." ] }, { "cell_type": "code", "execution_count": 1, "id": "abe920b9", "metadata": { "execution": { "iopub.execute_input": "2026-10-26T15:47:56.749199Z", "iopub.status.busy": "2026-10-26T15:47:56.748782Z", "iopub.status.idle": "2026-10-26T15:47:57.321186Z", "shell.execute_reply": "2026-10-26T15:47:57.320641Z" } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy.optimize import minimize\n", "from scipy.stats import chi2\n", "\n", "from cointf import ChebyshevKernel, Meioses, loglik_fast, simulate_meioses" ] }, { "cell_type": "markdown", "id": "e17f9ac1", "metadata": {}, "source": [ "## Data\n", "\n", "2,000 meioses of a 1.5 M chromosome, simulated with $\\nu = 4$ and $p = 0.05$." ] }, { "cell_type": "code", "execution_count": 2, "id": "c5ae6754", "metadata": { "execution": { "iopub.execute_input": "2026-10-26T15:47:57.322327Z", "iopub.status.busy": "2026-10-26T15:47:57.322202Z", "iopub.status.idle": "2026-10-26T15:47:57.331608Z", "shell.execute_reply": "2026-10-26T15:47:57.331127Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2000 meioses, 3077 crossovers, mean 1.54 per meiosis\n", "first three: [[], [0.519, 1.009, 1.451], [0.555]]\n" ] } ], "source": [ "nu_true, p_true, L = 4.0, 0.05, 1.5\n", "data = simulate_meioses(L=L, nu=nu_true, p=p_true, n=2000, seed=1)\n", "\n", "counts = np.array([len(d) for d in data])\n", "print(f\"{len(data)} meioses, {counts.sum()} crossovers, mean {counts.mean():.2f} per meiosis\")\n", "print(\"first three:\", [np.round(d, 3).tolist() for d in data[:3]])" ] }, { "cell_type": "markdown", "id": "ae57e458", "metadata": {}, "source": [ "## Likelihood\n", "\n", "`Meioses` validates the data and precomputes every distance once; each likelihood evaluation then\n", "reuses it. `ChebyshevKernel` is optional: it makes repeated evaluations several times faster and\n", "agrees with the exact kernel to about $10^{-12}$ in total log-likelihood here (checked below).\n", "\n", "The optimiser works on $\\theta = (\\log \\nu, \\operatorname{logit} p)$, so every trial value is\n", "in range." ] }, { "cell_type": "code", "execution_count": 3, "id": "b3aa8b33", "metadata": { "execution": { "iopub.execute_input": "2026-10-26T15:47:57.332398Z", "iopub.status.busy": "2026-10-26T15:47:57.332341Z", "iopub.status.idle": "2026-10-26T15:47:57.341737Z", "shell.execute_reply": "2026-10-26T15:47:57.341350Z" } }, "outputs": [], "source": [ "mz = Meioses(data, L)\n", "kernel = ChebyshevKernel()\n", "\n", "\n", "def to_params(theta):\n", " return np.exp(theta[0]), 1.0 / (1.0 + np.exp(-theta[1]))\n", "\n", "\n", "def negloglik(theta):\n", " nu, p = to_params(theta)\n", " return -loglik_fast(mz, None, nu=nu, p=p, kernel=kernel)" ] }, { "cell_type": "markdown", "id": "576164c1", "metadata": {}, "source": [ "## Fit\n", "\n", "Nelder–Mead needs no gradients and is robust for two parameters." ] }, { "cell_type": "code", "execution_count": 4, "id": "1fb81755", "metadata": { "execution": { "iopub.execute_input": "2026-10-26T15:47:57.342659Z", "iopub.status.busy": "2026-10-26T15:47:57.342605Z", "iopub.status.idle": "2026-10-26T15:47:57.672967Z", "shell.execute_reply": "2026-10-26T15:47:57.672570Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nu_hat = 3.773 p_hat = 0.0493 log-likelihood = -2698.526\n", "109 likelihood evaluations\n", "exact kernel at the MLE: -2698.526176 (difference 0.0e+00)\n" ] } ], "source": [ "fit = minimize(negloglik, x0=[np.log(2.0), -2.0], method=\"Nelder-Mead\",\n", " options=dict(xatol=1e-6, fatol=1e-8))\n", "nu_hat, p_hat = to_params(fit.x)\n", "ll_hat = -fit.fun\n", "print(f\"nu_hat = {nu_hat:.3f} p_hat = {p_hat:.4f} log-likelihood = {ll_hat:.3f}\")\n", "print(f\"{fit.nfev} likelihood evaluations\")\n", "\n", "exact = loglik_fast(mz, None, nu=nu_hat, p=p_hat)\n", "print(f\"exact kernel at the MLE: {exact:.6f} (difference {abs(exact - ll_hat):.1e})\")" ] }, { "cell_type": "markdown", "id": "40c99655", "metadata": {}, "source": [ "## Standard errors\n", "\n", "The observed information is the Hessian of the negative log-likelihood at the MLE, computed here by\n", "central finite differences on the $\\theta$ scale. Wald intervals on that scale are mapped back to\n", "$(\\nu, p)$, so they stay in range." ] }, { "cell_type": "code", "execution_count": 5, "id": "fa8a5c03", "metadata": { "execution": { "iopub.execute_input": "2026-10-26T15:47:57.673970Z", "iopub.status.busy": "2026-10-26T15:47:57.673918Z", "iopub.status.idle": "2026-10-26T15:47:57.697805Z", "shell.execute_reply": "2026-10-26T15:47:57.697460Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nu: 3.77 95% CI [3.30, 4.31] (truth 4.0)\n", "p: 0.049 95% CI [0.029, 0.081] (truth 0.05)\n", "correlation of log nu and logit p: 0.52\n" ] } ], "source": [ "def hessian(f, x, h=1e-3):\n", " x = np.asarray(x, dtype=float)\n", " k = x.size\n", " H = np.empty((k, k))\n", " E = np.eye(k) * h\n", " for i in range(k):\n", " for j in range(k):\n", " H[i, j] = (f(x + E[i] + E[j]) - f(x + E[i] - E[j])\n", " - f(x - E[i] + E[j]) + f(x - E[i] - E[j])) / (4 * h * h)\n", " return H\n", "\n", "\n", "cov = np.linalg.inv(hessian(negloglik, fit.x))\n", "se = np.sqrt(np.diag(cov))\n", "z = 1.96\n", "lo, hi = to_params(fit.x - z * se), to_params(fit.x + z * se)\n", "print(f\"nu: {nu_hat:.2f} 95% CI [{lo[0]:.2f}, {hi[0]:.2f}] (truth {nu_true})\")\n", "print(f\"p: {p_hat:.3f} 95% CI [{lo[1]:.3f}, {hi[1]:.3f}] (truth {p_true})\")\n", "print(f\"correlation of log nu and logit p: {cov[0, 1] / (se[0] * se[1]):.2f}\")" ] }, { "cell_type": "markdown", "id": "6b89dc8f", "metadata": {}, "source": [ "## Likelihood surface\n", "\n", "The log-likelihood relative to its maximum over a grid of $(\\nu, p)$. The outlined region is the\n", "95% joint confidence region (a drop of $\\chi^2_{2, 0.95}/2 \\approx 3.0$ log-likelihood units).\n", "The region is tilted because the estimates of $\\nu$ and $p$ are positively correlated: with more\n", "escape crossovers absorbing the close pairs, the remaining interference crossovers must be more\n", "regular, which means a larger $\\nu$." ] }, { "cell_type": "code", "execution_count": 6, "id": "0ccadfad", "metadata": { "execution": { "iopub.execute_input": "2026-10-26T15:47:57.698761Z", "iopub.status.busy": "2026-10-26T15:47:57.698709Z", "iopub.status.idle": "2026-10-26T15:47:59.258427Z", "shell.execute_reply": "2026-10-26T15:47:59.257942Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAG3CAYAAAC30lEJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAnUBJREFUeJztnQd4U3X3x79JN92lLZSWVfZGpgo4AEFAQMCBuLevKDj+LlSGOAAXrtc9eAUVEFHBCeJgL9l7Uwp00Ukn7f0/5xcSmrZpb5Kb3JHzeZ60yc29Nzf7m3O+5xyTJEkSGIZhGIZhfBSz2gfAMAzDMAyjJiyGGIZhGIbxaVgMMQzDMAzj07AYYhiGYRjGp2ExxDAMwzCMT8NiiGEYhmEYn4bFEMMwDMMwPg2LIYZhGIZhfBoWQzqjuLgY5eXltS6jPpo1recqZWVlKC0tdWmZ2pw7dw4lJSXQKnR89Fwx2kULr2stHINe3l9afawYbcNiSMYbSytv9szMTMTGxuKDDz6oddnhw4fFsi+//FKR273jjjvQu3fvOpeNGjUKAwYMgJaYNGkSGjVqBK3x1VdfoVOnToiPj0diYqJmXmNyUKJpvTP78NS6cqnpte5t3DkGT4oDLb6/tPB8MfqDxVAdXH/99WjSpAm0islkQlBQEPz8/Lx6u4GBgeJ2Gechsfrggw/irrvuEmI2KytL84/l5s2b8dBDD6Fjx46IiYlBq1atxJfOgQMHZO+D1qX3U1xcnNhHv3798Mcff3h1XV98rd90003o06cPfAW9P1+MOrAY0jn169cXX6YPPPCAV2/3888/x7p167x6m0Zhw4YNIr0wYsQImM36eAsOHToUBw8exKeffioE3E8//YRTp04J4bF79+46tz99+jQGDx4sohQ7d+7EyZMnRRRxzJgx+Ouvv7yyrqvwa11f8PPFuII+Pol1hFyfTuVwPqVI6MvRVZzxB9G6jlIyFRUViobe5e5PrfWcTbXQ/aX7Lcc3UfXxoX1aU65paWm2qB49H9Z9Wv1DdHImrSHn/rr7mIwfPx4///yzSD9QFLJ169aYM2cOCgoK8M4779S5/WuvvYacnBx8+OGHaNCgAUJCQjB16lS0b98ezz77rFfWdfV5l+OPc7RtbWk7Z15PNSHn9ULL6bm3+ghre/978rWkxOMl9zjq+mySe/xyj0XO/ngeuvZhMaQAFJl5+OGH0axZM0RFRaFp06YipUC/oCuTmpoqQtbk5yGvyK233ors7GyRc6fcuyvU5BmqCfoAvPvuu0XK79dff7Ut//PPP3H11VeLCFN0dDS6d++OefPmuZWXP3HiBEaOHCmOKyEhARMnTqzxw4lup1evXiK1Qbc9cOBAcTyurvf222+jTZs2Yp0ePXpg+fLlcCZac80114jngo6bHof33nvP7rjp/tJzVpUnn3yyWirV6p86cuSIiABRCoee+549e2LKlCliHTpPfqFHH31UXKYvb7pMJ3p90DYUkakpAkf3f9iwYbb1KELz/fffV1vHlee2JiZPnlwtimUVH8ePH69z+yVLlojnhLapDN2Hbdu2icfJ0+u6+rzX5o+jlOeQIUPEY9yyZUvb40upOzoG2mfjxo3xxhtvVLttZ15PNSHn9dK/f3/xPti3b59t3ar7lvs6cef9pcTjJfdYHX02yfkckfN6kHsccvfFaASJqZWRI0dK8fHxDq8vLCyUevbsKbVp00ZauXKlVFxcLK1Zs0Zq166d1L17d6mgoECsR/+7dOkide3aVfr333/FeitWrJDuvPNOKSYmRnriiSfqfCYyMjKk0NBQ6d1336112cGDB8Wyzz//3LbOwIEDpeTkZGnz5s229RYsWCCFh4eL2z516pS4L//73/+k6Oho6cMPP7StN27cOHHclalp2bBhw8R9puu2bNkilZSUSIsWLZIiIiKk1157zW7d2bNni2OcMWOGlJ2dLZ0+fVp64IEHxPH88ssvTq83c+ZMcTv0OOTl5UkpKSniOIYPHy4e39rIysqSGjVqJN13331SWlqadO7cOWn//v3S008/La1bt862Ht3fG2+8sdr2jz76aLXXiPWxGDNmjHjM8/PzxWNB0DHSfTpx4oTDY6qoqJCOHTsm3XPPPeLY6P5Y+frrr6WwsDDpoYceko4ePSqVlpZKW7duFfeXjt2Z59Ydli1bJu4H3f/aoOeN1hs/fny16+g46bqlS5d6dF13nvfaXus33XSTtHPnTqmoqEiaNWuWeMyXL18ujRo1StqxY4dY/sYbb4hjofd7ZZx5PdV0DHJfL/Qa7NGjR43byX2duPP+UurxcufzSs7niNzXg5zjkLsvRjuwGHJTDH388cfiTfbbb7/ZLac3cmWR8sEHH4jLa9eutVtv3rx5YnllMUQigj4UrKeysjKXxdCePXukjh07Sr1797b7gCQx1rRpU3H/qvLYY49JjRs3Fl+yzoohul36UKsMfRh36tTJdpmEQYMGDaQbbrjBbr3y8nKxz4suusip9ejDmZ6ju+66y249Wp6UlFTnh/U///wjjpvEbG04K4Zonxs2bKi2vhwxRB+e9Nzn5ORIUVFRtueXPnjpuRkyZIjDbZ15bl0lNzdXiPvY2Fjp0KFDta57+PBhcX8nTZpU7brff/9dXEdfJp5c153n3dFrnQQpfbFXfs6aN28uxcXFSdu2bbN7vbZs2bLa61NJMeTo9VKbGJL7OnH3/aXE4+XO55XczxE5rwe5xyH3tcVoB06Tucnff/8tKheuuuoqu+VXXnklwsPDsXLlSnF51apViIyMxMUXX2y3XtXtCApBU1jVeqI0kyvQbVMomPwdy5YtQ1JSku26TZs2iRQbhabJe0B5dmuu/ZJLLsGZM2ewZ88ep2+T0mJUcVSZdu3aISUlxZZbp9smrwnddmUoDUPL9u/fL3w1ctejSqezZ88KI21l6PHv27dvncdMj09oaCgef/xxLFiwABkZGVACCqFTKkwuZAamKrMWLVqIdCuF15OTk8XzQukFgh4Tem6oesoRnnpurdB+KL1Dhuo333xTHKMcavJNWJeRf8ob6yr5vDds2BAdOnSwXSYvFT13tM/OnTvbvV5p+bFjx6Akcl4vtSH3deLu+0uJx8ud17TczxE5rwe5x+GpzxTGc7AYchMyb5JgqelDl3L49OYgyBtE61WFllX1YlhLQ60nf39/l46Nqmno+G688UbxwVUZ65uT/Ankt6APKjqRmLn//vvF7dK2zlLVu0GEhYWJD4zCwkJx2bpf8jnU9JgR9LjJXY8e28rLalqvrmP+8ccfxf3/z3/+g+bNmwtvARmD5RjbHRko6bGUCxngyV+wdu1a0YOIPnDpvpEfLTg42Ga0pcsEHasjPPXcEvR43H777aJ0fdasWbjllltkiUKiptu0LiMfhyfX9cTzXtNrnb4AHb0HcnNzoZQhV+7rpTbkvk7cfX8p8Xi585qW+zki5/Ug9zjcfW0x3se1b1nGBn0g0xuEPsCqihr6tWH9JUTr0S+sqli3rfrrQwnIYPnvv//innvuEb+MqK+NFasw++9//yuMvUpR2y9xK9YvMWtVVdVfu9Z1nF2vpl9f6enpso6bDJc//PCDMJrT80RfMM8884wQcE899ZRYh3595+fnV9uWSsxrIiAgAHLZvn079u7di48//tjO/En3s3KHausHOpnxHeGp55Zep/RaItPyjBkzRFRCDhEREcK0S7/Aq2L9NU8NKD25rjvPu7OvdTnvAVdeT668XmpD7uvkn3/+cfv95e7j5c5rWu7niJzXgzPH4c5ri/E+HBlyE0qHUYiUyo4r89tvv4nQMlVzEFdccYX44LN+sFihfi2egsTZu+++iwkTJogTpTSsUPqGvlgXLlwIb0NVFfRFVrX6iX4xLV26VKTZ6NjkrkdpRVqPvqQrQ7/Qqj7edUG/7C699FLxuFHlzfr1623XUQqCvmQr/7Kj21i9ejXcxdo0s2rzzKpdxOm+0n3++uuvHe7LmeeWBI6c1gyUdiLx8+233+Lll18W1ZLO7PPaa68VIr+yiKP1rNVgVEHk6XVdfd49hTuvJ7mvF2v0paYKJrmvEyXfX67izueV3M8ROa8HV45DjdcW4zwshmRi7dFR+USKf+zYsejWrZsQG7/88ov49fT777+LvixdunSxlc7SrwiKElHIlHxGtB79aqAyWE93S33xxRdF75Xnn39elEgTdJv05qQS0XvvvVf80qQPN/pw/uKLL6r5A5SkXr164njIx0THRKXZ9Ov+zjvvFGX5r7zyitPrUWsC+rCjiAU13tu1a5e4X/QhVBeLFi0SbQfosaBfj/Tc0GNAv9CtYpagyBr9EqbSeFpvx44dotmlnNuoC/JVUanxq6++KkrCaf/0q58ie+Q1q5xCnT17tvilSekqOgZKkdBravjw4UKAOPPcUrqLfu3SPmvjsccew9y5c8Xr57777rN7H1T9oq1pn9Q+gFIH9NxZPRr/93//J7wt1ufR0+u6+rx7CndeT3JfLwR92dN7h0q9KSph7TMk93Xi7vtLCdz5vJL7OSLn9SD3ONR+bTEuoLaDW+tcf/31olqippO14oMqa5555hlRMUVVC1S99dRTT4nKjsqkp6dL999/v9SsWTNxopJgWodKVp977rk6jyUzM1Pc7vvvv1/rMqqyoWVffvml3faffPKJqP6hyjUqxSWoiuOOO+4QrQESEhJEm4AJEyaIKjQrdD1Vo1WmpmVUHjtgwIBqx01l9XQ81jYDVhYvXixK/hMTE0UVBm1ftdrOmfU+++wzUb5L9+OKK66QVq1aJSqNqMS1Nqhab/78+aJMmNoPUIUL3d7cuXOrrfvFF1+IChTa5+DBg0VJ+5NPPimOS85jQdBzRY9Hamqq3fIDBw6Iihd6bbRq1Up65JFHxGNGx0OVKpWh8lxal66j6pZrrrlG+vPPP+3WkfPcfvvtt6Lqhcr1a6ugcfQeoFOfPn3s1ne0z+PHj4tWEnTMVIU0dOhQafXq1TXepqfWdeV5d+a1Tp8Xl19+ebXlY8eOrfY4OfN6qukY5L5eqJqKtqfl9evXr1apJud14s77S8nHy9XPKzmfI858DtR1HM7si9EGJvrjiohilIHC+9TIbObMmSKaxDDehCIrZLTfuHGjy0Z9b+yTYRjGk3CaTGW+++478f+yyy5T+1AYH4TStM8995yiosUT+2QYhvEkHBnyItOnTxclliR8yPhIJmvKxVOvoZqMjwzDMAzD+IAYon4YW7duQW5ODpo0bYrWrdvUuQ0d8vbt20SKqXfvi21lkZUh09r+/fuE6GjVqrVT/TA8BR0vVeKsWbNGHB81SaPp2tSYi4yxDMMwDKM1MjMzhQGdzPeXX365+FFvNFQVQySAJk+2TJZu2qw5tm3dIsTNg+MfdrjNhg3r8b85XyAsLFQM95v2wovo2PFCLxG6O6+/NgtHjh5By5atxJNH+x07dhxGXjvKK/eLYRiGYYzAqlWrRIPPiy66SFQqUrU0DZylajwjoWpSf+7c/8HPzx+vzJglShaPHDmMJ594HD179UbPnr1q3IbWm/Ts8wgKDMR9991d7XoSQ3369MXj//ekrZnXyn/+xltvvYneF1+Mhg3ldwVmGIZhGF9FkiTccccdYvwPtW4g3n77bVHsQyNJaurqrVdUM1BTY7S1a9fgyv79bX12mjdPRps2bbFm9SqH23Xp0lWkl2prNHjJpX3supp27txFPKm1de1lGIZhGOYC1ND00KFDdt3mqX8SfZ/SuBEj4a9mDrKoqAhJifYdYhOTknDo4EFFb2vjxg1CJDVr5jjPaR24Z4WebGooR7ON5LbXZxiGYRhPUVOTU1eg77eq32tB52dhVmb37t22Jp+VO5pTZ3frdUZBNTFUXFRke2ArExYahqIiy0BPJaBuo3PmfC78QjUZra18t+hbLFjwje0yCaMlS39BSlq2aOWudfKKah7+l1/seGBjgYNtrGQV1f2mSy+SNwfpZJ7zb+CUXNff9Om5lteX0mTkyLu/jDrERQWrcrvxkSEev43GkcoXWTSKUG6f8SHKPPb1Q9w/prAQ977awoPlzxW0khQd6HEhFBIRC5SddXtfNAiX5lVWZsqUKaJTd2VohBQJJBr+W5no6Gjk5eXBSKgmhgKDLC8cig5VprCoEIEKjaeg1ufTX5iCHj17Ydy42qdrjx5zHYaPGGm7TE+0M1PHtSiE9AwLIcZVsaqGICLx7Q1BxPguIiJUdhZBF/0H8HPjO7K8BAVb3kdKSordD/2gGr53Q0JCRBES3Xblimf6fqQxJ0ZCNTEUGxsnmrKlpdtPEk5LO42EBMeeILmcPn0Kk59/Fh06dMRDD02oNlG+pgnjlaeMVx6eqFf0GhXSmhDiaJC+MKogoveF0tEhem8qFR2izwOlokNa+Ox0JTrkFfyCYPJ3XQxZy8dJCNWV9WjVqpX4f/ToUbRu3dqWNaGZbtbrjIJqBmoSQhd1645VK/8R+UsiKzMTu3butKsk27FjO1audG4y8unTp88LoQ54eMIj1SY7M8aEhRCjtoD1VHrWl5DzQ6wu6vqxx8jj4osvFj36aEizlcWLF4uMDlWTGQlVS+tvvfU2THrmabz80nS0at0a//z9N9q0bYvLLrvctg6VxdOE4X79LOMqSJFSw8XCs5a86fr164QviFQqNVekkN6Uyc+JyE7LVq3x22+/2vbVqVMnNG7cBEbCEykyvUaFlIYjQvrGiBEirUeHGGMRGBiId955B7feeisyMjJEn6H3338fzzzzDJKTk2EkVBVDiYlJmD37bfz991/IyckRvh0SQpUjOZ06dUaDBg1tlwsLzyI19YQ4f/WQoaJEny43aNBALKMoU4+ePcX5U6dO2t1eixYt4Cu4kyJTCy2lx1gIGQO1BJEn8YQgUgojpcoYCzfeeKMYJk5zNCnY8O2334oRUkZD9XEcWoUMYqSCtV5NpkYVmdaiQr4uhLKzlU/NREcbzwzsbVHkSf+QlivLtFJV5u2KMk9Xk1m/k4J6POKeZ+hcCUo2zUZubq6mv9u8DY+V1jFqpcjkwEJIP8JH7u3oXSB5O0qkt3SZ0aJD9KPPHUGkaRM1ozgshgyIJ1NkcqNCeosIaTUa5C3x48qx6FEcGTFtxjCMjqvJGO2hVlRIbbQkhEhwVD5pGb0cp5rPtyery5QuNtDa+1apzyOGkQOLIR9KkXFUSLtCSI+iQs/CyCiCSKt4IoLMMJ6ExZAPpci0FhVSOz2mthDSm4CQi17ukxEEkdGjQ+6i1cpZRnuwZ0iHqDF+Qyu/9IwghPQgFJS8n1r2FrGHSB6hgX5oGB6EAD8zysorcDq/BGdLyz387Fh+pCkxq8xV2ETtO7AY8hG88QvJG1EhPQshXxFBehNF3hJEnqouU7qyrHITxvr1AtGlUQQSI4PtppxTR5bU3GJsO5mHrMJSxavKTh4/jIK8PDRIbIz6IfYzIjMz0pB+6iQio6OR2LiZWJZy9LAYuZTYxHK5KseOHEJ+Xm615aGhYWje0jJmgvFtWAwZCLVTZHqICqkhhHxVBOlNFHkDPQ10bRIVgn7J9eFnviCCrJAwSooKQUJEMFYezsLxHGVf4x/OeB4b/lmGASNuwCtvfmR33cuTJmDVit9w9cgb8MIbluveePEZMWh0xnv/q3F/M6Y8hS2b1qNxFbHUvnNXTJv1jqLHzugTFkM6Q6spMk9HhfQohFgE6UcU6T1dpnR0iCJClYXQqZwyfLcpB8cyS9E0NhCje0QhISpAXE/r/bo33WGEyFWS23bEqmVLUZCfh7BwS3PAzPTT2LDqL7Rs28HpqPjFfS7D7I/nKXqMjHFgA7WPp8j0EBXSmxDSi4FYbbT2OHnrNaKH6rI28WE2IfTDvzkY8tpB/PePTPy0LU/8p8u0nKD1OjeKUNxr2LRFGzRr1Q6LFy+wLVvy7Tz06nsFYuMujGhiGCVgMWSQqJCrKTI9RIX0hJa+3PWClkSR2hWG7qDUeysy2B8t6tezRYSmfHcKZVW80nSZltP1RFJksDBZK83g0ePw++KvbT6lJQvnYsT1t7q0r9zcHOzc9q/dKa3K/ErGd2ExZHD0Xlqql6iQlr7Q9YpWHj9vvF60HB0ir5DVLE2psapCyAotX7zZEh2i9RuEuz4vyxGXDxmFYwf24ND+Pdi09h8xqLtv/6td2teendsx7amJdqflv/xY53ae+qHJaAv2DPkwSgxk9WRUSA9CSCtf4EaBHk8teIn06iFSwjsU5H/hNzJ5hGrjaKXrA/3MileV1QsNQ7/BI/Djwi9xJjMDQ0eNhb+/a19b7BliaoPFkAHQwy8XNdJjLIT0iVYM1p4WRFqtLCs5V2E7T2bp2mhW6frS8gvbKcmg0eMwfcIdKCkqxJdL/nFpH+fKJSgft2KMBKfJDIyeU2TuRoVYCOkfX4i6eSJd5u4PDyqTJ38OQVVjAQ6sQLR8VPcocZ7WT8svgSfo2O1itO3cHVdePQLNWnBPIMYzcGTIR0vqvZki83ZUyJNCyBe+oLWE2lEivabL3CG3+BwOZRWiZWyoKJ+fNjqhmomahBAtp+uJE7nFHu1IPe2/8+rsRJ2Xm4Pd27fYLavcmNFqoK5MYGAgWrfr6IEjZvQGiyGdo4cUmTd/LbMQMiZqeol8MV226kg2msXUg7/ZhJHdotArOVSYpckjRKkxighZhVB5hYTtJ/Nq3Z+zvqFGTZMREhrq8PomzVsgPCLSdrlxs2Rs3bgWM55/1G69PldchfsffVZcv2PzOmGarkx8gwS8N+dC6T7ju5gkazyUsSMvLw+RkZFISctGRETNPTS0EBlyJIbc6S2kZlRIi0KIo0HaQU0fkScFkSfEkLtG6tZxoRjRoYEQRI4gISS3A7WrJmor7s4oCwtx/bd/eLBF+NVGUnSgV76Tgq6aCVOA64+lVFaMkmVPITc3VxPfbVqBPUMMUwsshLSFUZ8PLZba7884i7mbU3Ew86zNQ2SFLqfkFInO00qP4tBig1gjRuAZezhNpgOcjQq5A0eFjP/Fq3fUSpnpzT+kRJk9Taf/dvtptIoNEX2EqHyeqsbSvDS1XkkoWu5OdIgxNhwZMiB6HL+htfQYCyFto9bz40lPmhajQ1ZI+BzOKsTe9ALx3xUh5M5oDobxNCyGGKfQWgUZCyHfhQWr93CmuSrD6BEWQzpFrRSZVjtNKwV/weoLNZ4vPUWHjD4LkGGUgsWQwfoLeTJF5olfh1pKj7EQ0if8vPkOWk3zM/qH3WSMqlEhLQghPX+ZZp8569T60TGOe7foGW+bqj1pptZi3yHrj6FGEe4Zst2ZU8YwnoTFkA5Ro8xTS1EhXxRCzooeV/ajd6FkJEGktaoyhjE6nCbTMFpKkRkRrQohEixVT2rcrh7R6nPqLFrw0R3fbT+6whEHt2/2+LEwjKdhMcSokiJTOyqktS9NLYoQLR6T1p5bTw8EVgpn09d/f/U+/vf0nViz6LNa1/v9q4/xxoRbsPTzd6EHKv9gvOO6Ifj6i49UPR5GO7AY8gG0aJx2BSMKIb2IDb0JIxZE7kWEVn7zgTi/Ys5bdoKo8mcBCaHvP3pDnP95znuyI0Ryfnw9d/+N+G7O+1CKCXeMwbxP7QVbbk4OiorkfxZwF2pjw2JIZ7gyi8wIUSEjCSG9CYuq6OXYtfBc6zFV1qR9N/S//cJA06qCqKoQIq697zG07NxdsWMoyM9DSXGRYj/uCvJzUeyE8GF8DxZDBvELeQqtRIX0/uWodwFUE3q4P956zvWSLpPLpWPuciiIahJCg8bdq9htz3zqARzcvQ3zP3kbN/ZtK075udl46s5R+PrDN/DGcxMx8rJOmP7UeLH+uKGX4vcli+z28dSDt+LTd2aJ888/ei/27NyKOR+8iat6JOPSDk2Qk31GXJeXm42pT03ANZd1w8j+vTDnI32k+xjl4Woyg+MN47ReokJqoXXBoNT903s1mhZRssze2aoyEkRWIWT9v/a7L1CUn+sxIUQ89PyrSD1yCN0uvQKjbn9ALAuLiEJBXg7m/vdV3PfEC3j48UkIj4gS1+Xl5qC0pLhaZKmoqFCcf+qF13Hs8AFc3K8/xt01HqHBfoiMihbXffb+W5g0/VXc+9Dj2LF1M5548E60bN0Ofa4YoOh9YrQPiyEd4e0UmdJRIV9LjxldBOlFFHmr5F4vpfbOUFUQeVoIEaFh4TD7+yMoJASR0fXtj2fAUIy85V7UD5Ev6sLCI+BP+wsOQVRMfbthrdeMugFjb7tHnE9s3BTzv/wU61f/zWLIB+E0mQbRSorMKO3+vS2E9JA+8rX7r3f/kNqCKCQ80m5ZaESkR4RQXSS36ajo/lq0amt3OSYm1pZCY3wLFkMGRku9hdSKCnnzS1CLIkBNtPZYeOO14ImUrpJGald+vJBXqHJEiDiblyu8Q67iauFGQKC8iJAkSbLWM5nNLm/LGAsWQzpHzRSZJ6NCehJCLIL4sTEqJISsKTKicoSITNTuCKLa8PPzQ0VFhax1Q8MjcLYg327Z6dSUKvvzhyTJ2x+jHKdOncKmTZuQk5MDrcNiSCfouceFGiXC3hRCjH4eJ71Gh7QghKi67PF5/9hVmXlKEMXGJ+D4of2y1u3YpTt++WGBKJ8/d+6cqBo7cfyI3TpxDRJw5MA+xY+TcQwJoD59+qBnz57466+/oHVYDBnUL+ROikzvUSFvfOFxNMi1x0wL6NE/5O1UWU1CyGqmrlp27wlBNPqO/2Dfjn8xqmczW2m9o8+3Bx57DoFBQRjUsyWu7t0KKccOo+NFPe3Wv/meh7Bz2yZc3ikRl1QqrWc8x7333otrr71WNw+xSeIEaY3k5eUhMjISKWnZiIiIUF0MOVtJVpsY8maKzNteIW8JIcY91K4480Z1mdKVZUpOsq+txJ46UNMojpqEUGV2Lp1j12/osbfnOtV4Uc70+uLCsygpKUZEVAzO5ucJzxBVhRFVK8pKS0rE9SaTSZTWUwVZcEg9u3WKaH/FxUhsFI/8vFwEBgYhOOTC4yrSbSYTQkPDajye8OAAh8eaFB3ole+koKtmwhTg+mtLKitGybKnkJub67Hvtg8++ABffvklfvzxR8TGxmLx4sWaF0YcGfLBknp30XJUyNOwEDLG46jH6JA3O1D3G/tArUKIoGoyKq8nht4+XtEO1FaC64WK8noSOGERkTYhVBMUHaL1rOX0VYUQEVIvVJTX03oRkVF2Qsha1u9ICDHy2LlzJ6ZMmYK5c+cK75de4D5DPpYi82ZvIWejQlpOj6n95W1E6DFVM0Lkrf5DWmzAWBeXj/sPmnftLYRRbZAgSu7YzSNCiIHHI02VCQoKEqeq7NixAyUlJQ73ExgYiM6dO4vzNOtt7NixePXVV9G8eXNdGKetsBhiDNFXiIWQPlFbEHkSvTdhrEsI0Q+nRhGBLgsh+mEmJ1XG2BMcHQ1TQPWol1ykskKQtGncuLHd8ilTpmDq1KnV1n/22Wdx8uRJh/tr0KABfvrpJ3F+8uTJqFevHtq3by+qyAoKCsTyQ4cOYd++fWjTpo1mn04WQzrFEykyvUaFWAjpGzUFkd6iQ0rh7GgOxnikpKTYeYaCaogKEeT7kYvZbBYtER54wJJmLS8vF//ff/997NmzB5988gm0CoshjeCscdoTKTK9psc8AafFvP94GzFCpHR0yJupMsbYkBBS2kA9c+ZMu8uUJouOjsZrr73GBmrGOGgxReaJqBALIXVQ63FnMzXDMFxNpkPjtNZTZHpOj7EQUhcjCiItRj4Zxhv4+/uje/fuIjqkdThNpmG0lCLzZF8hV2AhZFyMmjLTWqqMfUOMpwkLCxNGaj3AkSFGl1EhpeGIkLZQ4/ngdJl+PjsYRmlYDDHw9agQCyFtYqTnRUvCn2GY6rAY0hmu+IW8VUXmjS8HFkK+hbcFEUeHGMY3Ud0ztOXff/Hzz0uRm5uDJk2a4oYbxyI+voHD9WmU2vbt2/D7b78iNTUVDz00AS1btXJ7v0YuqfdGmNuZqBALIcYZ2ENk3BJ7dxov0mde1flkDKPLyBAJlldeeRHt2rfH7XfchbNnz+LZSc/g7FlL18qamPvl/7D4u0Xo0LETUlKOo7ikWJH9GhVvGaf1hpFSML6AN58vT0WHtJgqM+r7m2F0FRmaP/8r9OnbD6NHXycut27dBnffdTt+++1X27KqUISHOmVmZWbi008+Umy/ekCtwaxGiwrpWQjlZ+c7vC48Otyrx8IwDGMUVIsM0UC3gwcPolu3C3NtAgIC0LlzF+zcucPhdo5ahru7Xy2hlRSZVn41+pIQIrFT28lT2+oBI0SHGIbRJqpFhs6cyRL+n6rNmOgypb+8vd+ysjJxslJYWAi94w3jtDeiQkYXQt4SKlVvR4+RJL37h5Qcz2EU3xDD+LQYKi+vsByAv/0h+AcE4Nz54W7e3O93i77FggXf2C5XFkZa6Tztiyj1C11rQkgLkZrKx6AnYeQtQeSrQ1wZxhdRTQyFh1s+fPPz7b8UCvLzEREe4fX9jh5zHYaPGGm7nJeXh4SEBGjdL6RmiszTUSGjCSEtCKC6jk1PoojRBvRZ0iiCq7oYfaOaZ4jSVjExMThwYL/d8n379iG5RQuv75d8RfXq1bM76cEvpKfeQmqgBSGkJ8+OXo7VW8+rJ7xDaqeLteoNZBifLa2/atBg/LF8GU6fPi0u//XnCpw6dRIDB15lW+ebb77CzBkvK75fX8bbxmk1o0JqohdhYbRj9yW81fWdYYyOqqX1Y8Zcj4z0dEycMF6kt0pKSvHg+IfRvHmybZ0zWVk4deqU7fKGDevx1by5KK+w+H/ee/dtBAUF4+qrh+DqIUNl71dPKJ0i85UPYbWiQkYSEVpOn+ndTM0wjHYwSVR6pTIFBQXCoxMXFyfSVVWrw0jMWP071DgxKzOr2j4iIyMRGRUle791QdvRPlPSshER4bqHSYnO066IodrSZN70C6kVFVJDCBlJBDlCi6LIG4JIaSO1UhVlhBIVZY0j3fP8uOMZcrUDNeFqB+qwEOfjAOHBtX+HJEV71jdl/U6KvOEjmAJct3FIZYXIXXAfcnNzPfLdpldUH8dBhIWFiVNNxMTUt7scGhomTu7uV4t4wy/kbeO0LwghXxBBeogU+WqJPcMwBhFDjPOomSLTslGUhZD3RJFWBBGnyxjG9ygtLcXRI0eQk5sjegtWpmPHjrKDJlZYDGkcb47g0HtUyJtCyJeiQY7wJUGk5Z5D3HyR8TWOHzuG6dOnCRtNcHD1COv0F19BcjKLIZ/GnRSZkcuHlYSFkDYFEcP4AmFR4TAHuu4Zqij1Qy70zYKF89G5c2fcfc99irXBUbW03lepyTytlF/Ik/h6VIjLzR0/LlpACz2lnMHIPxgYxpPkZGej/4CBivYDZDGkQzzhF1Kyt5CzH/J6EUKMbz8+Ruh95QhfarzoSiUZoy1atWpdrbGyu/CrwkB+Ia2kyLwJCyHtoIVKMzZTM4zxGT3mOkybOhlZmZlontwCfn5+dtd369bdNppLLiyGGFX7Cmn517YvRDs8AfuIGG/gao8hRv8cOnQQKSnHxamm9jlNmzRlMWR0v5AeSuqNEBViIeT+46dWhMiT0SEtV5UxjK/w09Il6N9/AO66+16nGyo7gj1DGkVLKTItRoVYCGkfFpR1wyZqhnGesnNl6NvvMsWEEMFiyMfxJeOkHPgL3BiPp94qy3xhVqCeqWsUB+NdkpNbYM+e3Yrukz1DjCIYISrEQoiR9frjVBnDqEpSYhI++eQjMZA9uUUL+JntDdS9el/s9Nw1FkMG9wt5K0XmDCyEfAs2VDMMoyTr169DREQktm3bJk5VadW6NYshI6CkX8hIKTKOCOkXNQQRl9kzjDF5ZtJziu+TPUOMV43TWiul59QYP9ZqwyZqhlEfTpMZGCM1WvREVIiFkPHxVHSIfUMMox5//rkC2dlnHF5/5ZUDEB0d7dQ+WQzp3C+kdorMG1EhFkLGgf1DjJrwKA5jsGP7Npw4kWK3LDs7G2fOnEFCQiP06tWbxZCv+YU8hZFLdTkipP7jz5PuGYZxlQkTH61x+c8//4Tt27YiKamx0/tkz5BB8UaKTI9RIRZC2sCbz4Ov9RxilIV7DOmHoUOH4ciRwygoKHB6WxZDOkCrKTK9wUKIURKtFQMwjK9TUVEhTnl5uU5vy54hlfxCWkZOikxvUSEWQr6dLuMye4YxDrt27awW/SktKcG6dWsBmNCgQUOn98liSCVqMk8r5RcyQhUZCyHfgP1DF35cxEUFu/0jJj7SN4bI8sR63+brr+fhyOHDdsuCgoLQrFlz0YPIz8++I7UcWAz5GHWlyJQ2TqudSuCIEMMwjPfYvXs3nn32WWzYsAGlpaXo0KEDJk+ejP79+yt2Gy+++AqUhj1DGk+RKe0X0kOTOF83vBZnZ8GXYMHqu8SHuBcNY7RFeXk5rrrqKphMJqxZswa7du1Cz549MXToUBw9ehRahiNDjMfgqJDrwsdZQRQcXd+p9X0RT/iGuPmisXsMcSWZc5w4cQInT57EhAkT0LRpU7Fs0qRJeO2118QMsWbNmkEpKOp09MgR5OTmQJIku+s6duyI0NAwp/bHYsiH/EJKpMj0EhXSYrTBkxGfyvvWozBi7xDD6J/GjRuLSNAXX3yBbt26CR/P+++/j4SEBPTt21ex2zl+7BimT5+GM2eyEBxcPbo4/cVXkJzMYsjwVWRqpsg8GRUyqhDydtpL78KIYRjPk5eXZ3c5KChInNzBbDbjhx9+wJAhQxAZGSnSZQ0bNsTSpUtRv75yn0ULFs5H586dcfc996FevXqK7JMjQxrGkV/IE2ghKmQkIaQV34/1OPQgirwRHeISe0bvREXXg1+Q6+ne8hIg9XwUpzJTpkzB1KlTq63fpk0bHDt2zOH+kpKScPDgQXG+qKgIAwYMQLt27bB48WKEhIRg9uzZGDx4MDZu3KhYmiwnOxs3jbtZMSFEsBgyEO6kyIwQFdKCENKKCNKzKGIYxvOkpKQgIiLCdjnIQVRo+/bt1Tw5laHoj5Wff/4Ze/bswd9//424uDixbMaMGZg7dy4+++wzvPDCC4oce6tWrXHgwH506NARSsFiyEB+IUY9tCqC9CaK2DvEaKXHkNHN0ySEKoshRziTOrMKI39/e2lBfX8qiyZ3GT3mOkybOhlZmZlontyiWl+hbt26IzzcuSgziyEPoSe/kNopMj1HhfQigmo6bq0KIr3BFWXah6fVe4d+/fohJiYGjz76KN58801hbqY0WWpqKoYNG6bY7Rw6dBApKcfFKSysulG6aZOmLIZ81S+k5xSZu7AQcg1fFUTsG2IYzxAXFydSZU899RSaNGmCc+fOCc/RwoUL0atXL8Vu56elS9C//wDcdfe9CAhQJoLHkSFG91Ehb6PXaJBe0macKmMY/dK7d2/89ddf4jx5jZRMj1kpO1eGvv0uU0wIEdyBWmW/kNooOX7D2aiQHtNjRhJCvnC/GO3C3aeNj8kDQohITm6BPXt2K7pPjgypjDPmaVf8Qt5MkXkbbwohXxALWkqb+WJ0SIlhrWrQKMI7hmZvYnTztN5JSkzCJ598hIz0dCS3aAE/s72Bulfvi2WZwyvDYsjgfiFvpcjUiAp5C18QQloURAyjRCWZN8zTEV64DeYC69evQ0REpBjxQaeqtGrdmsWQkSvJlEbpCfVGiwr5kgjyNUHEJmqG0S/PTHpO8X2yZ8jAaDVF5m5UiIUQfEIIaqGJJsMwvgHH9jTWbNFRisyRX0iPKTJ3YCHkPXwhQsQwVWG/kDZZufIf5ObkoF+/y7B9x3Zx3hG0TmRUlFP7ZzHko6iVItO6V0gLEREtwYKIYWqH/ULeYeOG9Th+/Dg6d+lqO+8IWofFEKMIvhgVYiGkPXyxqswX8FZZPXeeNg6PPf5EjeeVgiNDBkWLfiF3okKGFULnSoCSXHEyleZZzpeXAqENIIU1AurFASZ1rX0cHdLHSI74SNdvv3Gk9svjvTWTjNE+Bw8cQMtWrRxeFxEZifj4eKf2yWJIB8NZlfYLKZUiU2P0hifwuBCiic/5KTCd2W8RPefFD0ryYCovqX1TcwAQlgApPBEIawQpLBEIJYFk31fD0xhREHFFmfsYrccQ+4X0wUcffSDGfdx9z70ICaknlpWXl2Px4kVYuGA+Zr36utP7dEkMlZWV4cjhw8jKyhSX68fGio6QVSfVMp7tL6TH8Rtaiwp5VAiR2EnbAlPavzAVZjhcTfIPBoIigcAISPTf7A9TwSmg4BRMFaVA3nGY8i7kxyWzPxCaACmiCaSkPkCwc0ZBvQkiTpUxWoX9Qurw4PiHMfvN1/F/jz+KiRMfQ1R0FN5+yzIQ9oknn0bTps2c3qdT6uXwoUNYuvRHrFu3FiUlJbZW2zR/hKbTXnzJpRg2bDiSk5Phq2ihx5DRUmS6EkIV54CsPTCf3gycOQATJFuER6rfTqS/SPhIQREWAUT//YKq7UZsJVUARZkw5Z8ECk7ClJ9q+U/RJIo05adAOrkeUuLFkJpcAQRYfiExDMN+ISPTrFkzEf2ZM+dzPPfcMwgMDETbtu3wxptvITo62qV9yhZD777zFjZsWI9LLu2DRx99HC1bthJ5OSIvNxcHDh7A5k0bMXXKc+jd+2KMf2iCSwfE2OPrKTJPRYUUFUKUBiORcnozTOnbYDp34bGRIppCatgdUlxHgCJAzkBeoXrxkOrFAw26VhJIZ2AqSIXp5AaYco/AdGIVpFObIDW5HFLiJYCf51IXRkyXMd6FZ5IxSnDu3DkUFhbCz88PFRUVqFevHgIDXR+jIlsMxcbF4YMPPxE3WJXomBj06tVbnG6/4y78+OP3Lh8QY7wUmatRIV0Ioay9MB/5Haazp22LJEp3NewGqUE3oF4sFEUIpFhIdIrrDJzZD/OR38Ttm478Bil1LaTkwZDiu9KURHgCFkSMEc3TzvqFOEWmHvv27cXs2W8gJDgEM2e9Dj8/M96a/QYee3QiHp7wCDp27OQ5MTR27DhZ65FYkruuL+Bus0Utpsj0bJxWTAiVFcJ0cCnM6VvFRcnkDym2vYgCIbqFdyrASOzUb4OKmFaWiNSRZTCV5MC0dyGk9O2oaH2tJRVnANg3xDCMlc8+/QS9evbCLbfejoAAi4h9+ZVZ+Oabr/DCtCl4/Y3ZaNy4CZyBHc8+hBqNFrUUFVJMCGXugnn/DzCVFUAiV1BS3/OeHZVKq01mSA0uEqk4U8pqmI79AdOZfTBvnA2pxVBIDXsoHiXi6JDy6HFivVbh/kLG5uEJE5GU1NhuGRVw3XLLbbjoom4ICHA+kuiSGCoqKsKSH3/A3r17UFBQ/Utr1qtvuLJbRuURHHJSZHqNCikihEoLYDq4BOaMHeKiVC8OFW2uAyLs35SqQSbtpleICJV53yJhsDbtXwwpYycq2o8D/KsbtRnG1R5DRiqr5xSZvqgqhCrToUNHl/bpkhj66MP3sXv3bvTt2w+hYaFwl8zMDOTk5KBRo8QaPUmubEMVbmlpaTh7tgCx9WOdbs2t10oyLaGVqJDbQkiSYMrYIYSQqewsJJghNe4HqVl/IUA0R2g8Ki66H6YTq2E6ugym7AMw75qLik63i7J9peDoEOMNuNmiPVFRIfAPdj0Kfa64AkYyUWdmZqC0xD5A0DAhQVSYOYNLn4ybN28S+bmkpCS4Q2lpqTBBbfl3M+Li4pGRkY5bb70dQ4dd49Y2R48exRuvvyqEUExMfaSmnkDXi7ph4sRHERSk7q9jtfxCdaXItNZbSDOU5sN84EeYMneJi1JoA0s0iJogahnTecEW1RzmbZ/AlHMIpr0LILUbq3pHay3BjRe9C1eSMUpADRY/+eQj/LF8mThflZmzXhMV7x4XQ0FBwS7X8ldmwYJvcPDAfvz3vx+KirT169dh1sxXRJvt1q3buLwNRa7i4uLw5uy3RdldRkaGcJn/+svPGHntKKhpntYzzqTIdB8VomgQmZIpGnSuCBKJiyZXWLxBCkZXPE54Eio63ALzjjkwU7osYAmkliMU8xB5MzrEJmpGC34hriJTn59+Woq9e/bg5Zdn4v3338O4cbcIz9DXX88VxunmzZ3vdejST8RevXoJReYuf674AwMHDhKihqD+RE2aNBXL3dkmLy9PqEISQgQJo/r164vlesHbfiEj4roQqrB4g/YusAihsARUdHsQUrOB+hJCVqJbQmp7vTB7m0+uh+nYCrWPiNHQXDJfh0dw6I8d27fhxrE3iSAIfc+HR4SjS9euePa5KaIpdGlp7WOOasKlT/YzZ87g119/wZo1q9CwYYKtE7WViY88JmMfWcLzk9yihd1yunNHjhx2a5uxN92EOV98gbj4eCGEtm3bipLSEgy+ekitI0boZIWaOfkKSqfItBAVcl0IlcO07zuY07aIixVNB5yPBnl3FpjSSPGdz7cD+BHmY3+gIjAMUqPeah8Wo1O8ZZ5mvxBTE7m5ucILTISGhiI/z/LdER4ejpiYGKSnpzs9ksMlMVQvNBRXXNkf7lCQX2A7+MrQ5fyCAre26dSxM1q1aoV5c78UESEyUo8aPVoII0d8t+hbkYKzUlkY6cE87Y5fSA6eriLThBCqOAfznm9gytxtMUm3HSNK1o0Cje2oKCuA+dgKmA78CCkwHIht7/Z+2UjN+IpfiFNkWkECzsdgSPT8/fef6Nyli5iZevr0aURHWzJHHhdDDz88Ee7i5+9nM0RXhlzhjga+ytmGqshefPEFIXw++vhT0ZCJ3OZPPfkEKsorcN31N9S479FjrsPwESNtlymllpCQAL2Yp7WCFmeQyaK8BOad82DKOQjJ5IeK9jcpIhS0htR0ACpKC2A+tUGU31eQEdwgjRkZ75XVG8UvxCkyfdKqVWtbFTkVT01+/lmMu+kGMZZjxMhrERER4fQ+VTNAxMbGwWw2i5Rb1VRYXGycy9vQdYcPH8JN4262daak7Xr27IWNGzc4FEO0rnV9a8meWuZpb/qFtNBbSPWoUFkRzDu+sAw+NQeiouMtwmdjSEwmSC2HQ6LZZvmpMO/7DhWd7nDbUM3RIYZhvMW99z1gOx8fH4/Zb72Do0eOICw8rNYeRLXhco3t/v378OqsGZg48SFMnDBenD9wYL/s7anEvU2btkKgWCkuLsb27dvQuUtX27JTp07h8OHDsrcJCwsTgikry/5LMTMr0yW1qAc8nSLTCy4JISqd3/axRQj5h6Ciy93GFUJWzH6oaHM9JLO/6EFkOnXh/aR1PDWvjtEeWvMLcYpMuwQHB6Ntu3YuCyGXI0Nr16zGG2+8Jtpe9+t7mcjd7d+3D5OeeQqPPfZ/YrK9HCh6Q3NE5s79H9q0boOff/4JkZGRuOqqQbZ1Fn/3Lfbv3y+Un5xtSDANHHgVvpr3JSoqytGgQUNs27oVW7f8i0nPPu/K3WU8mCJT6svNJSFUnA3z9s9gKsoS/pmKTncCYQ3hE4TGQ2o+GKZDP8F06GdIUS2UHyrLGBJXzNPe8Atxisy3KC0tFdGgnNwcYY+pTMeOHREaGuZ5MTR//je4+577cHWV6iyqMKPr5Iohaps9bdqL+OXXn3HwwAFRIj9h4qMICblQJpqQ0MguZSVnm3vuvV/kFLdu3YJ1a9cK/xD1I2jdpubeRe4ixzxdk19IbbSQIlONwnSYt38OU0kupKAoS0QoxDv9crSClHgJpKw9MOUchnnft6joei9gcr1qjlNlDMN4g+PHjmH69GnCIkNRoapMf/EVJCd7QQydOnUSl112ebXl/fpdhs8/+8SpfVFoi06OGDV6jNPbUN+B/gMGipOWqck8raRfyJspMl1Fhcgrs+NzmMoKLfPFOt/lmyZik1l00zZvegumvOMwpay0tBFgnCY6OkS1Ia2u9hjSsnmaU2RMbSxYOB+dO3cWQRm5I7w84hmKioqq0R9Ey+g6Rht+IV/AJSFE4ylICIUloqLrfb4phKwERwlDNWE6+gdQpMBAW4bxMpwi8y1ysrNFsEMpIeRyZGjQoKvx5huviRI2SkcRB/bvxw8/LMaIEdfC1zHKGA5PpshUMcJSSfmuuTCVl0CKbIaKjrcB/l7sfXLmpOV/TCNoCeqlJNHoETJTH/8LUpvq0Vi5cKrM2GjVL8T4Fq1atRbBF1cn1CsmhqgnT3BIML5fvBjzznwpltFA1BtuHIuhQx0PWWWM2XVard5CTkWFKsph3v21xSMUUh8VHW/1vBCyih+5y6viLdFkMqGi2QD4kRg6vQVSkyuBEOebljGMEeEqMu1BGmTa1MnIysxE8+QWttFbVrp1616tObNHxBCN3xg2bLg4WcdWKBmu0hOumqedabboyC+k15J6JaJCzqbHTId/gSn3CCS/QDG4FP4emuUkV+jI3Ze3BFFEE1FRJqbbn94MqflV3rldhnHTL8QpMt/j0KGDSEk5Lk7UTqcqTZs09Y4YqoyviiBP4Mg8reUUmR46TpvStsCcukacr2h7PRDaQLsCyNG+vSCKpIQeFjGUtgVSswHCYM0YEy2bp7UER4W0yU9Ll6B//wG46+577Zolu4NsMTRjxsvi/9NPT7KddwStw/jOYFZNR4XyU2Hav1icraD0T2wH6EIEqSCKpPrtIfkFwVSSA+QeBaKSPXZbjP6m1bNfiNEKZefK0LffZYoJIafEEA08rek8oz3ztJZTZF6l7CzMu+bBVHEOUkwbS7RDbyLIm6LILwBSXCeYTm+yRIdcFENsovZeWb3R4BQZI4fk5BbYs2e3Ogbqe++933a+V6/e6FJpZEZltm3bCl/BG5PqXfELaTVF5tWokGQ1TOdYDNPtblAm7aOmEPKCn4gqy0BiKGOnpeTej9MpjG/CKTLXoGGpv/76Kz777DNkZmbir7/+qrbOoUOH8Mknn2DHjh2IjY3FqFGjMHLkhUHpdZGUmIRPPvkIGenpSG7RAn5mewN1r94XOz1+yyXPEI3DWPTdD05f54vINU97yy+kZorMm+kx0+HfRGdlMXhVCcO0VkSQpwVRZFNIwdEwFWfDlLkbUoOaf/Qw+sVbfiGtldTzhHrvMHDgQDEWq3Hjxli0aFG169evX4/bb78dd9xxBx544AExbuvWW2/FI488ghdeeEHWbaxfvw4REZHYtm2bOFWlVevW3hFDjsjPz7cbi8EYFy0bp01pW2E+sUqcr2h7nXuGaS2KIE8KIpNZRIdMx1ZYUmUshjzefdqofiGtpcicgaNCrrNw4UJhpZk7dy4+/vjjatd36NABu3btsiuHpwr15557DlOmTKlWJl8Tz0x6Dkrj1KuJJtPXdJ6gQWlU5ta2reMxGUbCGykyNfxCnmq06G6KTHZUqODkBcN048uBuI7GFUIeEkQiVXZsBZB9ECjJA4Kc+4VFsG+IYXyT+nV4imsqhafJFWVlZWIOqRwx5AmcEkOVw05VQ1B+fv5CCF3ZXyGTqg5R2jyttF9Iryky2UKorPC8YboMUnRr93rl6EUIeUIQhdSHFNEUprxjMKVvg9S4H7QECevwaOd6iDgiOiZUkf0w9nCKTPvk5eXZXQ4KChInb1NSUoK33noLgwYNUuX2XRJD9z/woPgfHhGBceNu8dQxGRpnmi16Gy0bp+tEqoB5z3zhdZGCY1DR7kbXDdN6E0IeEERSXEeLGMo+pDkxZCS4kswCp8jkExsZhAA3/Fhl57Ma5OmpzJQpUzB16tRq648dOxanT592uL+4uDiRGnMFyijdfffdSEtLw9KlS6EmLiVdWQi5bp7WerNFvUaFTKlrxWwtyRyAig43AwEhviWEFIZmtwnyjgmhqaUGjEpFhXwRV8zT3vALeRI2TtdMSkqKXYYnyEFUhozNxcWOv6NcjeaQELrvvvuwbNky/Pnnn0hKSoLuxNDGjRuwds1qTJj4qN3yt996E5f26YsePXrCyBjVL6TbqFBJHkxHl4uzUouhQFiC7wohpaJDYQ3F6BIaaouzaa4/pozuzdNaS5GxcVoZSAjJqbi6+OKLoTQkhKiS7Mcff8SKFSvQvn17p/dRVFSIoKBgmM1m4TdaufIfBAYG4JJL+rjkO3Lp595X8+bi2lHVJ1vTsq+/mufKLhmV/UKenFDv8ajQoZ8sk+jDG0NK6Om7QkjJ+2LyE/PKxNncY+7vj2EUSpEx+ufBBx/E999/L4QQVZc5y9EjRzBt6hQhqoh33nkLn3/2CT784AN89NEH3osMnTp1UjRKqslFfvJkKnwRLXSe1hpeiQqdOQBzxg5IMKGi9UjX0jlGEkIKIkzUVFFGqbJE5X8dGgG1yupdgVNktcPl9MowY8YM0XSRfEDEFVdcIf7Pnj0bXbt2Fd6gDz74AMnJyRg/frzdtvPmzUNiYmKdt7Fk6Y8YOnSYiAAVFBRg44b1eOfd/4rmiw899B/cfvudTs9NdUkMNWyYgE0bN+Cyyy130godUHx8vCu79Ak8aZ7WYorM41Gh8jKYD1gafEqJlwJhLqSHjCqEFEiXSZFNbZEhy+8vRkmMaJ7Wc4qMUYZrr722xtRas2YWH2Lv3r2FR6gmYmJiZN1GWtppxMVbqoV37dyBli1bITY2TlyOi4tHVlYm6tWzRLbl4tIra/iIkfjww/eRejIV7dt3AH1S7t69C0uW/IC777nPlV36LN4wT+sxRSYHU8rfMBWfgRQY4fzcMaOKICUFUURjSDBbBrcW5wDBUUoeHWNAv5CWUmRsnFaHtm3bilNt1WfWaJGrxNaPFaO/2rVrj7//+RudO3cRyyltlpWVVWPmyiNiaMCAgcKwtOjbhfh24QKxLCamPm67/U5xnZGRa56WW0mmlXlkWkqRyYoKFWbCdPzvC6Zpfyd+kfqCEFICvyCLcbogVZTZS06KIW68qG+0VkXGxmnGytBh12DqlOfx09IlIlV299332sZ0dO7cGSEhzqXICJdjjldfPUSccnJyYKKRRlH8q9HIaGr8hiTBfPBHmKRySNGtxKR1xjPRIUqVmQpSATJRx1t+fTH6wxvzyJxNkbFxmnGV1q3b4J1338fxY8fEoNbIyEixPDQ0FHfedbdL+3Q7AUtttBl1qc0vpEbXaU9HhUwZ24WxVzL5o6LVCBpsI/8GOCrkHKENLY85pSNhLLj7tHFxJkXGxml9QgVbVUd/dOrU2eX9uSyGqMZ/9+7dyMzMQHl5hd115PL2JeRUkmm187Tu/ELnimA69LM4KzW5XIyOkI2vCiE3okNSkOUXF0pylT0mH68kc9U87S2/kC+lyBj98eefK5Cdfcbh9VdeOQDR0dFO7dOlV9ihQwfx8kvTxXlKk1nd2xUVFUhIaORzYshV9NJ52pkUmcejQkeWwVSaDykk1iKG5OKrQshdgs+LITJQM7rE11JkHBUyPju2b8OJEyl2y7Kzs3HmzBmhQXr16u0dMTTni88xaNDVuHHsTRgzeiQ++PBjZJ85g9mz30Dbdr4xtb42fN087THyT8B0cr04K9JjZv616HGCLGlw0Yn6XLFzRnWGUQCOCjFVqTr9wsrPP/+E7du2IinJfu6axzpQHzlyGMOuGS7Om0wmUVkWHRODB/4zHn+u+MOVXTIq+IW0lCKrMypEg1j3fw8TJFTEdwWiW8rfOUeFXMcvEJL/+dQMldgzqmGUFBlHhRhPQVkp0ifUiNErYqiwsBBhYWHiPLm4KUVG0JyTvLw8GBWtziQzQoqsLkwn18FUcBKSXzCkFkPkb8hCyP3H4Xx0CMXsG9Jbs0UtpsgYxlOQVYdOeXnOf1a5nWdo174D5s79UpTZ//3Xn2jSxNK11ldw1TytF7+QJqJCNIj1yDJxVkoeDATKnFrOQkg539DZU6L5olEqytytJNPTGA5fSZGxV8h32LVrZ7XoT2lJCdatW0s/ndGggaUK1uNi6Kabbradv+222/HmG69j2tTJaJiQgEceecyVXTLsF/LsIFbGZaiiTDQv4Ioyw+NsikxLxmnGd/j663k4cviw3bKgoCA0a9Ycz0x6zqWp9S6Joeuuv8F2Pj6+AV6ZMQvl5eUuHQCjXb+QN1JkdUaFco+5NoiVo0LKldhb02Ql6qfAw6NlRgUNhit+IW+kyDwJR4UYR7z44itQGpc8Q+PHP1BtGQsh9yvJmOqYj/8l/ksNu8sfxMpCSFnM5xvYVdSdEmZ8bzirXDgqxGgZlyJDBfn5OHu2AKGhFhM1w3gkKlRwCqYz+0RUSHZPIRZCynO+w7dJqjCMZ8gdjOoX8nSKzBk4KsTUxfHjx7Hl383Iyc0Rw+Irc801w1HfyWGtLomhSy7pg19+/tkuXcZ4xjztif5CWkmR1YVtEGtcR3mdplkIeQZralIyhhTS2xgOI6TIOCrEKMnmzZsw45WXkJiYhJj6MTBZXI02Skqd/950SQyRElu27DesWrUSSY0bI8DffjcTfcRELaeSzN1mi3qZR6Z4VKgoC6aMHeKsU52mGeWx+bTsx+4wxkmRcVSI0RO///YrbrhhLK6/4UbF9umSGKLJsFdc2V+xg2D0i8eiQin/iAaLUkxreV4hjgppLjIUHO3E3DiGYRiZUKPn9h06Qklki6H169ehd++LxfmHH56o6EEYBaOYp51JkXkE6it0+l9xtqIxR4XU53wIWuLIkLf9QlpMkXmynJ69Qowc2rRpi927dqJDhw7wuhiaNfMVLPruB3F+7I3X4Zv538KXMEr3aS2M4KgrRWY6sRomqRxSRFMgqnndO+SokGfL6zXiGdJzWb2RUmSMusRFhCCwnuuivDSgXNHjUYNGiYn44P33kJmZgeQWLeBntm/r06v3xWIihkfEUHh4uJgSSwPQKETFuIYS5mmt+IU8kiIrK4Lp1PlhrOwV0lQ1mRE8Q3ozT/uScZqjQoxcVq38B+HhEdi2bZs4VaVV69aeE0P9BwzEo49MEKKIuPuu2x2u++lnc5w6CF+qJFPCPK1n6owKnVwLU3kppNCGQEybunfIUSEvpsmMUU2mF7wxmFVLxmmGkQt1mVYa2WLottvuwKWX9sHp06fE+I3bbr8TvoySlWTeQumSesUpL4Updc2FCjJbRMIBLIS8w7nzr5sAY/bX8bRfSMspMk/CUSFGTzhVTdayZStx2r9vPy6//ArPHRWjeTyRIjOd3gxTWSGk4GhLbyFGG5RZBLIUYOwUk6+lyDgqxOiJlSv/QW5ODvr1uwzbd2wX5x1B60RGnR8j5MnS+rvuvseVzQyN0pVkWvcLKZ4iqyiHKWWlOCs17geY6phzx1Ehr4sh+NeDns3TevILeSNF5kk4KsQozcYN60XX6c5dutrOO4LW8YoYYjxvnlYDNVNkpoztMJXkiOiD1KC7asfB1EBZoeW/D0eGtJ4i01o5vdpEhPBXm9F47PEnajyv6qBWRjm8ZZ5WsqTelRRZrVEhqeLC6I2kPoDf+cGgjuCokFeh1KUgQL3IkC+hReO03qJCDOMsLJ9VriRjSNzsh6kwHZJfEKRGvfkhMYBniLtPGweOCjFapLS0FEePHBHjwaQqla4dO3Z0epA8iyGDVJK56xdSM0VmtkaFSAj51/GrmKNC3scAkSF3/EJGS5FxVIjRO8ePHcP06dNw5kwWgoOrv8+mv/gKkpO9IIbKy8uxdMmP+OOP5UhLO435CxaJ5Z9//imGDRuO+Ph4+FL3aS2Yp72F4imy3KMw5R2DZPKDlHipewfHeNZArZJnSM+dp33NOK027BXyDRYsnI/OnTvj7nvuQ7169dTzDP34w/f4/fdfMfLaUTh37oIQaN48GQsXzlfkwIyGmuZpLYzgcIQ55R/xX2rYDQiqo2MoR4W8z7limCrOR0J92ECtVbRUTs9eIcZb5GRni0bQSgkhl8XQ8uW/Czf3gAED7ZZ36tRZlLwx+sKTKbJao0KF6TBl7YVE8+mT+ta+IxZC6lCYIf5JgeGAfxD0iF5SZBwVcg+OCvkOrVq1xoED+xXdp0tpsszMTDGjjDBV6hIcGBiI4mJtl48bsZLMW/2FlG60SANZBfXbAvXiFN03owym82KInx/twVGhC7AQ8i1Gj7kO06ZORlZmJpont4Cfn31fum7duttGh3lUDDVs2BAHDx5Ahw4d7cTQ6tWr0LhxE6f2lZ6ehhUr/hDdJJs0bYoBA64SosrdbchdvmnjBuzYsQNBQUG4/IorkZSUBK1XkjnyCxkuRVZ2FqbTW8TZCmqyWBscFVI/MuSEWFWyksyX/EKMPVxKzzji0KGDSEk5Lk5hYdWN0k2bNPWOGBoxchTefutN3Dh2nLi8Zcu/2LrlX/z66y946OGJTjnCn332aXTp0hWtW7fB77/9hr//+gsvvvQK/P39Xd6GfEyvzpqBo0ePYPDgIQgMChLHO/6hh9G0aTPoES2Yp5VMkZnSt8MknYMU1giIaOrV42J8JzLEKTL9eoXkwlEh3+OnpUvQv/8A3HX3vQgIUKb/lEuvTPIKVVRUYP43X4v/L06fhpiY+rj33vvFTBC5fPnlHDHr7P+eeEpcvuzyK/DA/ffgr7/+xMCBV7m8zZIff8Du3bswe/Y7qB8bK5YNHny1Sym8/KJzCA9wXFavdCWZVv1CiqfI0ixRIanBRbUPZOWokLoUnY8MhehTDBkVT3ecVhuOCjG1UXauDH37XaaYECJclulXXTVInHJzc0VKKsrJOSBlZWXYvn2bEFBWaB8dO3XGpk0baxRDcrf5fdlvQiRZhRBBD5qSD5xWKsm0Oo+sduN0Bkz5JyDBDCm+izcPi3GGinKgKEu1yJAeU2RGME5zVIjROsnJLbBnz25h1VEKt2KWJSUlyMvLFedDQkKEN8cZEzals+Kq9CSKj4sXd9LVbc6ePYv0tDThNl++7HccOnwIMdExQkUmJCQ4PB4SWnSyUlh4vtGcjsdwaNUvZI0KIaYVEOhcYyzGixRnwSRVQDIHAkGRuus8rUaKTO/GaU91m/bE2A1Oj2mXjIwMzJs3T3xnv/jiizWuQ1mlX375Bdu2bUPTpk1x/fXX1+kXtpKUmIRPPvkIGenpSG7RAn5mewN1r94XIyKijlYtSoghEkHz5n4peg1ZBQTdiUGDr8bNN98q6w6VlpaI/1W7R5Kosl7nyjbFRRYBsGD+N2jVujXatW0nSvAemfgQnnt+iij/r4nvFn2LBQu+sV2uLIz0ap7WJDSHLG3rhRRZbXCKTF0KMy9EhWpLZTIuo7WokDM44xWSC6fH9M8jjzyCb7/9Fi1btsTff/9doxg6c+YMrrnmGqSnp2PMmDFivbfeegsbNmyQdRvr169DRESkEFJ0qgp993tFDJEi27lzBx56eILw7xBUXTZv3lwhRv7z4EN17iMkpJ4tklOZgoJ81KsX6vI29UIt6zRKTMSjjz4uzl89ZCiKioowf/7XDsUQleoNHzHSdjkvL6/WSJIRzNOe8gvVmiLLO26ZTk9zyOq3c2q/jALENJK9qqngpPgvhRqro7yex28YPSokF44KaZfhw4dj5syZWLhwoRA5NfHggw8Ki83WrVtt1WC7du2SfRvPTHoOSuNS08U1q1fhif97Cn37XoaGDRPEic7TMiqvl0NsbKwQMCnHj9stP378OJo0aeLyNiSY4uLiq+2DSv6zshx/SZOfiLpZVj4xymNKt6h4KbZD7dPpOSqkOqbc8++zCOfaZWjBL+ROioypG44KMY4YMGBArZaZ06dPC6H09NNP25XFd+jQAWrikhiitFRV3w5By+g6WTdsNqNPn75YsWK5iNoQ+/fvEykt8vdYWfHHcnz11Vyntrn88iuwdcsWW6qLfEZU/t+yRUsYCa2apx0ilcOUsdNyNr7mCB2jEaQKIN8ihiQVxJDe8IZxWitRIS6lNwaU/ah8Kimp2Z6iNJs2bRJ+od69e2POnDl4+eWXRVrNVWuKUrgUv+zevQe+/XYBbrvtDlvnRxreSsu6de8hez8333Kr6CL52KMT0KRJU+zatRNDhg7DRRd1s62zd+8e7N+/H+PG3SJ7G0p5kUiaMGE8WrRoiSOHD4leQ3fceTfcwZWy+poqyfRknlY0RZZ9CKays5BoxlV0C/cPjvEcZ9NgKi8V6UyENtCdedqIxmkjw6Zp+SRGBCA41HW/VrG/JSLfuLFlioSVKVOmYOrUqdXWf/3115Gdne1wf5GRkXjiiSdk335OTo74f+utt6JFixYii/Pcc89h2rRpWLNmjdPNElUVQ7l5eVi+fBnWrV0jhrNSaf2RI4dFGqpnr954a/YbtnUnPvKYw/3QnZ456zXRE4geIDJfU0fpytAwth49ezm1DYXoJk+ZJgRRRkY6rrlmuKguq9qy25fx5DyyOlNkcR0BUy3PBafIVMeUd8xyJqIxYHIpgOwzKTKOCrnnFWLTtDqkpKTYmYyDHKS2aHnVoqWq1zuDNTU2aNAgTJ8+XZx/5plnhDB65513MGnSJOhGDIWGhuKKK/vbLevU2bV+MSRQHJmaibZt2zm9DUFjQtq0aStO3sITlWR67DxdI+VlMGVa2h9wbyEdkOt7KTKjRIWc7TbtCdg0rX1ICMmpuHrooboLopyhXTvLd/oll1xiF11q3749Dhw4ALVwSQw97MTIDcYz6M4vdGYfTOUlkKhfTW1fsBwV0kYl2fnIkCRzVIoRUmRGiQo5dSwe8grJhdNjvkebNm3QpUsXrFy5EkOHDhXLKA23c+dOjBgxwuF2f/65AtnZZ2TdxpVXDkB0dLRTx6V8zSOjOlr0C5ltKbLOXk+7ME5SkgdTcTYkmCxpMi+itxSZ1jBqVIjRD9999x3+/fdfIW4I8gMRd999N5o3by7Of/LJJxg8eDD27dsnPEM//vijiBjVFoXatXMHTpw4YbtM1hwyYlP1eFBwkGjASIVViYlJ6NXrYu+JoaKiQuzevRuZmRkoL6+wu27o0GGu7pYxol/oXDGQtU+c5RSZDsg7X1JPxml/z/fO8aUUGUeFLsBRIWMSEBAgPEY9evQQJytUDW6Flu/duxc//PAD8vPz8f777wsPEdlbHFF5CPy6dWtFY+Unnnza1g+wtLQUCxfOx6GDB5GUlOT0cbskhg4dOoiXX7IYn8jETMosKytTqLSEhEY+IYa8UUlmFL+QKWuvZUJ9SCwQVksjS06RaQJT9kHxX4psBj2hB+O0HqNCniilZyFk7KaLw4cPr3O9uLg43HPPPS7dxsp//sa4m2+xa4xMky+ooOree+5Cfn4ewsOd60DtUr5izhefY9Cgq/HpZ3PE5Q8+/BgfffQpOnbshD59+8KIVC2rVxIjjOGoLUVmytgh/ktxnXisg9aRJJjO7LecjWntVb+QWoNZOSqkDJweY7wFda8uKa4ebKAWP2VlpbbyfY+LIcrVDbvGovworEXNkqJjYvDAf8bjzxV/uLJLRmPmaWf9QrWmyKxfriSGHMFRIW1QmG4Zl2L2B6KS1T4azaLnqBCbphm907XrRZgz53PRTJmaKhNpaWl4++3ZovFzo0aJ3kmT0UR3a68AKomjFBmN5KAyPepk6Yu4U1bvTfO0t/1Cpqw9MEnlkGjYp8zmfYx6lWSmMxZvlxBCfuqbcT2dIvPFqJAn4PQY402uHTVaVJa98vKLwp7j7+8vgjLJyS3w9DPPutRT0G3bf7v2HTB37pe4+uoh+PuvP0VXaIaxYhu/EUuNFnnyudYxWY3uMW18IkXmChwVsofTY4y3IfFz730PYMx11+PY0WMoLi5Cg4YJolqtNhN2rft0ZaObbrrZdv62227Hm2+8LkZkNExIwCO1dJz2dbwxhkNTJfWVU2TxnCLTPPR85R5zSgxpAa2X03syKqQF07Rc2DTNKE1MTH1xUgKXxNB1199gOx8f3wCvzJgljEs87qL2SjJnMEIlmV2KrB6nyDRP9gGYUGF5vkJivHazejJOay0q5AzOeIXkwukxRi1Ic6xetRKHDh9CcRFFhhri8suuQP3YWJf2p1h3LCMLIVcGtHqqkkxPnadtKbLaqsjYOK0hv9B+VVJk7sBRIfVM05weY9SCehNNnfI8Tp06ieQWLRAcHILt27dh4YL5mDDhEVxyaR/viKGNGzdg7ZrVmDDxUbvlb7/1Ji7t0xc9evR0ZbeMDs3TslJkNJiV0TZShdMl9XpGi1EhLZimOT3G6IHFixchKCgQH3z4iW2+Gg2M/2npEnz44ftiYDz5ijxeWv/VvLm4dtSYastp2ddfzYOvoZVKMi2V1FNVkkiRUaNFTpFpn5wjMJXmQ6KO015stqgn47SnvUJaKKWXC6fHGDXZv28frrv+RrtBs2Scvmb4CAQFBeHkyZNO79MlMUShqdga8nL169fHyZOpruySMRrWCfWxHThFpocUWdqWCylN6jFk4BQZR4W8FxViGE9A3aYL8qv/cKeeQzSfjK73ihiinkKbNm6otnzjhvWIj493ZZeGx6fGcFScu5ByiW2v9tEwdVFeClPmLnFWiu/qtceLo0KuwVEhxtfp2asX5sz5QmiOkpIS0WsoNfUEXn9tljBQN2zY0DueoeEjRoq8XOrJVLRv3wGQgN27d2HJkh9w9z33ubJLxk3ztKb8QjmHYSovgRQYAYQ73wmU8b7RXTxfwTFeTZEZNSqklVJ6NU3TXEbPeJLBg4cgPS0dr702S0SDqICLqstatmyF//u/J13ap0tiaMCAgaLb46JvF+LbhQvEMqr1v+32O8V1vo7csnqtzSRTzC9kS5G1A0wuBR8Zb6bITm8S/6WGPWQ1xlQiRWbkqJDeUKt6LELh/TG+g9lsxu133IkRI0bi6LGjorSe+hw2a+blposEdZymEw1Eo5uOjIpydVeM0aqSsvZYztavJUVm9JL6syeB0EbKrecpCjNhyj0KCSZIDS+CHuCokDqmabnIjQoxjLvQTFQ6KYFLYqi0tBS7du3ERRd1Q1RUlEiRffrpx6Lp0Q03jEVAgHHeDAUlZQgLDHHYY8hIlWSKpMjOnrZUJdFcq6jm8EVMaethTt+IioaXQIrr7ni9jM0wn16LiviekBr09uox2o4h7V/LmZhWQFCkLozTRsFTpfScHmN8gXItNF1csOAbMaiVxBANbZ0542XhHVq3bq0QSnfeebdLB8N4xjztzeGsppzDljPkPZFRlaQlmjSMxmXdkhEZFozcgmL88+9hHD+d7dxOzp4UQogQQkdUaHV3KITEeukbUR7WWNkIkZwUGUXxTlvEUEVDx6JNadxJkRkpKuQMnhi7wT4hRq/ka6Xp4sp//sErM2aK89u2bkFiYiKeenqScHPTjDIWQ96dSVaXedqrfqGcI+K/RFPPdZIi694uCZPuHoCrL2kjctFWqELh17X78PKnf2DznhPydhbaSESEbEKnBkFUWQiJ22l4iTqpsjP7YSrNg+RfD6jfrs7VfSkqpKUGi2pXjykJ+4QYrTZddEkMFRTkIyTE8ito584duOgiywd9bGwcCgoKXNklY4QxHJJ0YdCnTqqSRl7eAXNeuAlBgdXfCiSMhvZphwE9W+H2yV/jh78t5ed1QcKHBFBNgqgmIVRbKs2TmE+sthxvw25ei+LpJSrkLJ6MChmpeoxxn4TwQISEuR4pLDIHGrrp4pIlP4imi02aNHFqny6V+jRt2hSLF3+HvXv3YNWqVbioWzex/PjxY2jSpCl8Gb1WkiniFyrKgulcISSTPxCmoinYiYhQZSF0KqcM7y3PwJPfpIr/dJmg62k9Wl8uQhBRxOc8JIDMuz/2jhCSkyIrOAVTziFIMENKvHCcjH5L6eXC6TFG7wR6oOmiSz8HqYT+tVdn4rtF3+KqqwaJ2n7ixx++Fz2IGN/ElGeJConeQjrwC1FqzCqEfvg3B1O+O4Wy8gvXf/xXJqaNTsDIblFivWfuGoDrnpgje/9VI0TUy0cLESFxLNaoUFwHyn9pvpxey1EhPaXHlO4yzf2EGDWbLlKGqnOXrqJoi/xDc7/8n3ebLrZt2w6ffPqFMEtXVmAkkmoa08GohzfN08g7Lv5JkU017xciszR5hAiKAFUVQgRdpuW9kkOREBWAIZe2QZOGUTh+Okf27ZDgkTL+tRNCkl+Q54SQnKhQaT5M6dssx5LkvNFQL+jdNK2H2WNyYJ8QY9imi1aqhqLi4uJgZDxdVu9KJZmmzNNWv1CEc7laNaCqMatZ+rtNOdWEkBVavnhzDh4cECfW79ctGfN+Pl+OLgPyCFUWQmJZeYlYrlZkyHRyvWWIbnhjQMZzpdeokNbw1fQYwxi26SIpsKVLfsQffyxHWtppzF+wSCz//PNPMWzYcJ5PZmAc+oXOFcFUmG45rwMxROXzVo5l1i5qj1a6PjJU/q/wqmZpighZhVFtZfcepaJMiCFxPBwV8lpUSC/pMR63wegJ1Zsukjdo+fLfMXrM9fjve+/YljdvnoyFC+dj/PiHFTk4I6BEWb0uKsnyUsQ/KaQ+EBgGrUN9hKw0ja39i6JZpetzz8p7LhxVjdn1F1JaEMlIkZnStsFUdhZSUJTFL1QHvhIV4vSYfNgnxGiB48ePY8u/m5GTmyPmo1bmmmuGO9180aVqMhJCjz3+RLU5ZJ06dRZTZBnjVpK5lSLTiF+IoIaK1EeIGN0jCgF+Na9Hy0d1t4yaofVX/nu+qWQt1FY+X1OVGa3vFSQJphOrLGepgszk4E7rHF8xTWt59hj7hBhPsnnzJjz+2ET8+ecKHDt2VFSyVz6VlDr/HevSqz8zMxNJSY3F+cr5OfIQFRfrIIqhclm9t/Bq5+nz5mk9pMgI6ixNDRWpjxCZo6lqrKqJmoQQLafriV/W7KvbPE0dqOson6+pD1F5vQT3Gi/KMU5Tk8XCdDEqRUroWefqHBWqGU6PMYy6/P7br2L01/U33KjYPl0SQ1S2dvDgAXTo0NFODK1evQqNG+vjy5Bx3jzt0C9EzRYLUi1nyZSrE6izNDVUpLJ5Kp+nqjEyS5NHiFJjFBGyCqGS0nN45bM/5HWgju9Z52yyyoKI1vd4B2pJgvnYCsvZhF6Av7Y7FWs1KqSn6jH2CTFGpaysDO07dFR0ny6JoREjR+Htt97EjWPHictbtvyLrVv+xa+//oKHHp6o6AH6CmpVkilC8RmYzhVDorRLaDz0Ao3YoM7S1saLJHyoaqwqJIRoPbkjOWjoqpxZYySI3I4IyY0K5RyCKT8FktkfUuN+da6ux6iQN0rp9ZQeUxL2CTFaok2btti9ayc6dKjb9+hRMUReIfJPzP/ma/H/xenTEBNTH/feez/69btMsYNj9GGeNuVbokIIa6iLZouVoREbAx74QDRUpD5CVWeTUWqMIkKyZ5NZkStwvDST7EJUqCcQ6LpI8QZGME0brXpMDuwTYjzJhg3rkZebK86Hh4dj3rwvkZmZIQa1+pnt/Y+9el9sN6pDDi6/E6jzNJ1yc3PFgLSoKIvJ1Kh4useQrjkvhqSwRF2Yp6tCQoc6S1Mjxn7dmovyeaoaI7O0Mw0WVUFWVOgwTLlHReROanyZ5qNCRjBN+1p6jIUQ42n+WL4Mx46dn3IgBFEEtm3bJk5VadW6tefE0JEjh0XpfFUiIyOrLTt8+DCSk2uZWu4jeHpavVYwnfcLiTEcOoZM1fN+zobRuBAV6gEEVX+/agmOCtUMCyHG13lm0nMe3b9sMTTjlZfQqlVrDBp8tTBOU/vrqo0Yt2/bhmXLf8ehgwfw4UefeuJ4dY+3yurlVJIpZp7Ot0R9JJ2LId0hJyqUewymnMPno0KX17m6L0SF9JgeUwPuMM34ErLF0Ftvv4fvv/8Ob7z+qnByU5SIokLU6ygnOxtHjx5BUFAQBl89BBMmPAJfw9tl9ZoyT5eTedofqNdA7aNhHEWFGnYDgj2fytaDaVoLaN0nxIZpRmusXPkPcnNyhC95+47t4rwjaJ1IJ607ssVQcHAwxo4dh1GjxmD7tq3Yu3cvMrMyQYX17dq1x6jRY9ClS9dq88oYY2NvnjZmEz9de4WyD0Aymb0WFdI6WokKyb59NkwzjIAaOlPXaZpSbz3vCFrHY2LICkV/evbqLU6MemX1mqHgfIosrJEuzdOGRaqA+fAvFyrIQpSZ3+NLUSG102PsE2KYC9DUi5rOK4W+6qB9DF2U1Z9NuxAZYryDnBlkpzeLqB0Nh5Wa9td0VEirpmlPoJZPiCvHGG+SkZGBXbt2ITY2Fu3atavmMSbIbrNp0yaUlJSgR48eCAtTd6YliyHGPc6LIYn9QtqhrBCmw7+Js1KzgV7pK8SmaWVROiokBy6hZ9zl9OnTmDBhAv755x+0bdsWR44cQb169TB//nx07tzZtt7u3bsxbNgw0deNegZRyfxXX32FIUOG1OkZkoNHPUO+jFI9hgxXSVZeAlPJ+RdnKJunNRMVOroMpnOFkOrFQ2p0seGiQpweU88wzTC1kZWVhRtuuEGIHxrVde7cOXH5pptuEpEiK7fccgu6du2K7777Tqz33HPP4eabbxbiqaZ2PURdPiGve4YY9XsMaaaS7Gy6+CdR5CGgntpHwxD5J2E6uUGcrWg1ok5Tu5ql9FpMj/miT4grxxiloPEYlUdk+Pv748Ybb8TYsWNRWFgookQ7duzAli1b8O6779pmmz722GOYOXMmlixZIoRSTXjCJ1SZC7MHXIByfSkpx8WJzjO+BU1AF9SrZR4Zm6e9FxWiYawHf4QJEiriOgNRnm98arT0mB58Qjxqg1GCvLw8u1OJh77DV62iAe6NhRAirB2ju3TpYlsnJiYGTZo0qbGbtLdwKTJED9q8uV/i999/FSYogkrqqSHjzTff6nPl9d7uMaQ5v5COhrMaOj2WvhWmvOOQzIGQWjjOvVvh9Jg+y+jlwIZp4xIbEoxQN14rZ8st39kkUCozZcoUTJ06tdr6K1asEFEdR4SEhGDAgAE1Xrds2TK8//77+Oyzz2zLcnJyhEYIDbWPDNevX19cJ5f09HT8/NNSpKaewJ133Y1GjRKxedMmJCYlomHDBHhFDH3yyUfYuXMHHnp4Alq2bCWWHTx4APPmzUVxURH+8+BDruzWJ3FUVq+rSjI2T6vPuWKYDp0vpW96pVfGbng7KuQsvpwekwMbpn2blJQUu/ldQUFBNa43b948pKWd/6yvAaoYq0kMrVu3DmPGjMHTTz+N2267zbachBAFUchPRGk0KyS45AZSTp8+hf97/FG079ARhw8fsok1MmR//tmnLo3ucEkMrVm9CtOnvyymxVohJdYoIRGTJz/LYshXOJ8mk9g8rQHT9B8wlRVAComFlNQHvh4V0lN6TKs+IR6+amxICMkZZvrpp86P1lq/fj0GDx6MBx98EC+++KLddc2bNxfD3VNTU9G0aVOxrKKiAidPnhTXyWHx4u8wZMgw3HzLrXjyicdtyy/q1g1vvvmaGCDvyIitqGeIwmJx8dVTI7SMrmPUQ+lKMoeISrLcuj1DjOfJT4Upda04W9HyGsDsbyjTtB5nj7EQYnyVjRs3CiH0wAMPYMaMGdWu79Onjyinp0oyK3/88Qeys7Nx9dVXy7oNSo11vaibOG81YVuhKjLal1ciQ92798C33y7AbbfdYWumRINaaVm37j1c2aXh8VZZvdI4LKsvsrzYJP8QIIAFsHqm6XKY9y+GCRWoiOsExLSGr6fH9OQTkgOX0DN6Yd++fRg0aBDat2+Pfv36YenSpbbr+vfvL0zUdKJo0VNPPYXS0lIRwXnhhRdw5513omPHjrJuJyw0DJmZGeJ8ZS2UfeYM0tPSUL9+jHfEUG5eHpYvX4Z1a9eIga0U8jpy5LDoMUBjOt6a/YZt3YmPPCZrnxQmo3yfM8jZho6N8pMk2mrqgimnx1CYcT77laP4vPIOjlb7SHw7PXZiDUwFJyH5B0OiqFAdcHpMfz4hOXAJPaMFMjMzReSH+OCDD+yu6969u62ijBozUkps4cKFoiBr2rRpuOuuu2TfTr/LLsdX875EYmKiuFxRXoH9+/YJPzPNSA0Przv9p4gYIhf4FVfat/jv1PlCmZwzfLfoWyxd+iPy8/ORlNQYd919Dzp16qzYNmSm+umnJRg2bLhYz13kNFz0hR5DpuIzljMshtSj6AxMR5eLs1Ly0Do7TXN6zDW0bpiWWzkmB/YJMe5AQqhyNKg2hg8fLk6u3U5fnExNxbOTnhZG7GeffVoER9q0aYMHxz/s0j5dehc9/PBEKMFvv/2CRYsW4smnnkHr1m3www+L8fJL0zH7rXfRoEEDt7fZtGkjtm/fhkbn1aMn8Nmy+vORISm4lnAk9xjybE+hA9/DVFEGKSoZUsPuMFJ6TEtdpo0ghNgwzRiJsrIyXH/DjRg0aDD27d+HkuJiJCQ0QstWrcR1Z88WIDQ0zHtNF91lyY8/YsCAq0RYi4zXN954k3C3//7br25vQym7Dz/4Lx555DEEBuir75EuyuptaTLnWp4zCvYUyj4IyeSPilbX2ifODZAe01sZPZfQM4z3oO/2goICYZbu1au3SJtZhdBrr87EqVOnnN6namIoPz8Pp06dtGvdTa5w6huwb/9et7ahcNlbb72B4cNHopnMUj0t9RgyTGSI8QwleTAdtISipWb9gXqxHn+kjVI9ZiQhxCX0jK/i5+eHl196wa5rNgmhV1+dKXoiJSQk6EcMWTtNRlTpBRAZEelwMq3cbb5duEAYq4ePGCn7eOiBpMZNlU+MAySJDdRqpsf2fQfTuSJIYYmQkvrVuUt3o0JGEULO4CuVY+wRYvTIPffej+CQELw6a4bwDAkhNGuGqCSbNm260ykyTQxqraiQ7C5LUkWdIf/atiFH+c8/L8WMma/aRoVQRRlFi6iMz1GHSzJlL1jwje2ydVs94bUeQ+eKYCo/r8jZQO3d9NipDTBl77ekx9pe55VBrEZB6z4hFkIMI4+AgAA8+eQzmDrlebz91psoLi5GRkYGpr3wotPNFlUXQ9HRlvRKXq59FIg6R0ZHRbu8zZGjh8UD88jEC45yUo40THb58t8x76v5NZbYjx5znV0kiQbXuRJq84keQ+ebLUoBoYBf3b8+GYU4exqmQz9bHvvmgwAvdP42SlRID+kxOXBEiGEsBAcH49nnJuO5556B2WTG1GnTXRZCqoqhsLAwJDVuLGacXXKppS8BRW/o8pX9L8w5oWaOtJyUoJxtBg8eIk6VefyxR9ChQ8daS+tp/3SqLKAYB5Sdj0CRGGK8ExUqL4V59zeW6rHo1pCSLq1zl5we048QUqpyjFNjjFF56cUXcPTo0WrLS0qKYTKZxawyKySSmjVrpp802ahrR+PDD99Hx06dRX+A7xcvFlGdqyuJGXKN79+/H7Pfekf2NmriEz2Gys77qQIsDbQYL6THDi6BqTAdUmC4JT1mMmtSCLmCr0eEWAgxTN0MvnqIyNjIwWsdqJWCGjeSkJk390vk5uagSZOmmDxlGurHXqiO8fP3t4vYyNmmKgEB/vDzd777tCsNF32hrB4shrwrhNK2wnx6MyQautHuRiAwzGfL6PVimGYhxDDK0qNHT3gS1Q3UVw8ZKk6OuP/+/zi9TVVmzHwNnsBnGy6eK7zgGWI8S2EmTAe+tzzeTa8EopJ92ifkDHKjQp7oMC0HjggxjHZQXQz5KrruMWT1DPlzmsyjUaGKMpj3fA1TeSmkyOaQmtqPwKkJ9gkZp3KMPUIMc4HJk5/FkcNH8ML0F/H555+K846gdWhuqjOwGGKch9Nk3kmPHfoFpoJTkALqWdJjHvYJeRNf9gl5c94Y9xFijML119+Igvx8xMXF2847gtZxFhZDBsFrPYbsDNScJvOUEELGTphPrhNnK9reAARFGGoIqzP4mhCSS12VYyyEGCPRqdIw9rqGubuCqrPJGJ1yXgxRxILxAEVZoss0UdH4MiCmtccfZqP4hIwghJRIj7EQYhjn4MgQ4zwV56vq/PQ1AFcf/YRKYN41F6byYkjhjSE1u6rOXbJPyDOGaT03VWQYo3qG5MCeIQ2i1+7Tdc4mI+rwsDDOzh2rgHnPQpjOpln6CXW42ePjNnytn5DRS+g5IsQYlevr8AlVhj1DKuMLDRctWGfD1T5DjnHBMJ212zJ3rP3NXvEJuYqz6TEtCCFnkBMVYiHEMN7DEz6hynCajHFuLpldZIjFkGJCKHUNzKmrLQ8vdZiObOLxVyYbpmuGhRCjZWKCAxDmhhE/6Jzx0qwff/QBhl0zHI0aJbq8D85zyIS7T1eGI0OKVo5l7oHp4E/ibEXzwZDiO/usT0gPlWMcEWIYbXHgwH4UFp6vcnYRFkOM81gjQ5wmc18I5Z+Aec83MEFCRUJPSFQ95qNCyFMobZiWA3uEGEZfsBhSYRSHrrtPCzhNpogQKs6Gecf/LkyibzWiztSjkYWQmoZpLTZVZLM0w8gjMCgI5jqKTeqCPUOM83BkyH0hdK4I5h1zYCorgBSagIr2NwEmz1aOeRMWQo5hIcQwyvLii6+4vQ+ODDFuRIb4wXNJCNHMsZ1zYSpMhxQYgYpOtwH+QYbpMM1CyDEshBhGm3BkyEdGcSiK+Xw1Qnlp3cLgzEn4DHKEkFQO8+75MOUegeQXhIpOtwNBkR4/NF+qHOPUGMMYm0nPPIXDhw/VeJ2fnx+io2PQq3dv3HDDWAQHy/uMYTHEOE9gGFB8Bigt4EfPKSEkwbT/+wu9hDreCoQl1LmZUX1CLITsYY8Qw8hj9Jjr8Mbrr6JPn77o2KkzgoOCkJqail9//RkdO3ZCi5YtsXTpEuTm5uLhhyfK2ieLIR9BqSGtgoAw8c9UWmArsvdp5Aqhw7/AfHozJJgtHqGo5Do3YyHEESGCO0szzAW2b9uGUaPG4Pobbqy0FOh32eV44v8exfiHJqBnz154+KEHce+998uKDrFnyIMYchQHfa9TZIgo48iQLCFEwvH4nzCfWGV5/NqMBmLb1bkNCyHlu0tz1RjD6J+jx46gbbvqn6FxcXEICwtHWtppxMc3QHx8PDIzM2Xtk8WQDkZxaA6rGJKTJpMpFgwthI7+AfPR5eJ8RfJQSA271bmNWpVjWuolpHRTRRZCDGMMwkLDsG7tmmrLDx86hPT0NITWC0V5eTny8vKFIJIDp8l0gLbmklVKk1FZOHwUuamxY3/AfGyFuFjR/GpIjfvWuZmalWPO4qnKMRZCnBpjGEeMGj0Gk59/FkeOHBEzy4KEZ+gE1qxZjasGDUZkVBR+++0XDBo8GIGB8nqIsRjSyCgOPSEFno8elHq5ik1vQujocpiP/ykuViQPgdS4X52bcQm98YVQXUZpgj1CDOOYVq1a483Zb+OHH77Htm1bUFRcjAYNGmLCxEdx8cWXiHUGDx4CZ2AxpBEyclzvaK1WZMgnPUOyhdAymI//JS5WtBgKKUnbESEt9RKSi9wSejmwEGIYfdGwYQLuv/8/iu2PPUNeRv+jOMgzdP4LtyS/UjdqH/ANyRJCFTAd+rmSEBrGQkjjvYRYCDGMfqEBrVlZWZDkfBfVAkeGvDyXzBAER0Ey+8NUUQoUZQL14mBo5Io5qRymfd/BnLZFXKxoORxSoiVka5SIkLOwELLAqTGGUZa1a1Zj7rwvcfrUKXGZvEEDB16Fm8bdgnr16jm9PxZDjMMv6OLsrJqvNPsD4UlA7lGY8o5DkiOG9NqNWq4QKi+zTJ/P2iP6CEltx0BqcJHhhJCRmipyRIhh9MmmTRsxe/YbGHbNcFx0UTfRRyjl+HEsXrwIaenpmDTpOaf3yWKIcQkpoglMuUeB3GNAw+7GfBTlCqFzxTDv/NIyYsPsj4p2N3mljxDBQkjfZmk2SjOM8yxfvgw33XQzrh012s5U3a17D/zngXuRnZ2N6Ohop/bJniHGJaSIpuI/RYZkoxfvEB2n3GMtLYB52yeVZo3dyULIB4QQiSAWQgzjmNJSeVXYJSUlyMnJQVlZGeSSl5sjRm5UJSoqCjH16yMvLxfOwmKIcY2IxuIfTV5HWaFxHkVnBFvBKZj//S9MBSchBYSiosu9QFTzOjfjiJD+hVBdcESI8UWKi4vx9ttvo127doiIiBDi5M477xQGZ0dC6JJLLhFRnJ9++kn27SQmJuGfvy1FKpXZu2cPzmRlyW60WBlOkzGuERgGKaQ+TEVZQF4KUL+NvO206h1yNmqVuQvmPQtgqigTj0NFx9uBerF1bsZCyAILIYYxHps3b8bhw4exePFitG7dGseOHcOoUaNwyy234Jdffqm2/hNPPIE2bdpgyxZL0YlcRl47Gs88/QSOHj2Kzl26IDgoGCknUrB+3VqMHn0dQkLYQM14OVVGYsiUdwySXDGkRZwRQtRDiOaMnR+vIUW3tHiEAuqOcrAQssBCiGGMSZ8+fcTJSvPmzYXgue2221BUVISQkAufk0uWLMGyZcvwxx9/4JtvvnHqdpKSkjD7rXfx3XffYveuXWLfDRMS8MSTT6NHj54uHTtHhhjXiWwKpP1rqShzZjutRIecjQaVl8K0bxHMGTvExYrESyG1GAKY/OrclIWQBRZCDONbHD9+HGFhYXaT41NTU3Hfffdh6dKlLpXBE/Xr1xcT6ZWCxZAC+NSQ1ioVZQJKk50rAvxD9CGIXDFyF2bCvHseTGfTIJn8ILUaASlB3i8QFkLaFUI8XoPxVfLy8uwuBwUFiVNN61VUVDjcj9lsFv6gmqCU2axZszB+/HiYTCaxjAaojhs3DhMnTkT37t2FeVoLsBjyECm5PiCQ6sVDqhcvTNSmU5tlDSFVFVer2TJ2wLzvO5jKSyAFhKGi/U2yjNIECyF1hBCXzjNGRU4lY62cs2zbuLGlCMbKlClTMHXq1Gqr9+3bV0R3HJGYmIhdu3ZVW56WloYhQ4agd+/emDZtmm05iSNKa1FkiIRQbq6l8uvs2bMoKCgQUaSamDz5WRw5fETWXXxh+oto3jwZzsBiyEegfjT52fnKNV4kTCZISX1g2r8YppNrISVdCpjMrokTT0WJ3CnnryiH6fCvMKeuFhelyGaoaDcWCKr5V1BVWAhZYCHEMNojJSXFLqITVENUiNi+fbvT+05PT0f//v3RpEkTfP/99wgIuNBva+/evdi/fz+Sky1ixTpGg8TR3LlzazRaE9dffyMK8uV9h8XFcTUZ42Wk+K6QDv8GU3E2kLUHiO3gnmhRShS529OoOBvmPfNtfZQqkvpBSh4kyx9EsBCywEKIYbQJCSFH6S13yMjIEEIoISEBP/74o51XiJgzZ47dZYoOUWn9vHnzcO211zrcb6dOneFJODJUB2eKyxAapvzEbcPgFwCpUS+Yjv8F84nVqHBVDLkbLVKqoSMNWj25HiYSeBWlkPyCUdH2OiC2vddEEMGdpdVJjXEPIYZxHeonNGDAACGyvvzyS9FHiE4ELSN/kVZhMcS4jdSoN6SUfyzjOfJTgfBEZR5Vb3esLsqyeINyj9haBwghFCJP4LAQci4ipLU5YyyEGMY9Vq9ebfMXUePFylAvISq1rwoZqyMjI+1SaWrAYohxzzdEBEVCiusEU/o2mFLXQGp7vb4eVYoGpa6F6cjvliaK5gBIyYMhNbpYtgeKhZCxhRB3lGaYuhkxYoTT1WEkhLRQUcZiiFEEKfFSgMRQ+nZIyVcDga5NUvc6hRkw71tk8wZJUcmoaD0aCImRvQsWQvqePM9CiGEYFkOMMkQ0tkyyzzsuPDdSs4HafmSpUuzEapiOLodJOieGrJKIkxJ6iSo5ubAQ0q4Q4qnzDMPIhcUQo2h0SIihE6shNejmVHTFq+QcgfnAD5Yhs2KkRitUtB4FBEd5XQQZwSzNQohhGL3DYshJ0ouKPfNMGAApriOk1KZiVpl599eo6DCOVAM0AfWyyD4Ic8pKmHIOWhYF1IOUPMQi3FSIBhEshC7AzRQZhlELFkOMMiZqwmRGRdvrYd78LkwFqTBvfAtS86sgJV7iXDNGpdNhGdthSlkF09lTYpEEM6SEHpCaDwICnJuLo3ZazJWIkDPRILUiQuwPYhhGTVgM+RCudKF2mpAYVHR7EGbqSp17BKZDP0FK346KNqOB0AbwGueKYTq1UVS3mUos7d4lc6BFBCX1cTpipYVoEMFCyDFslGYYxlVYDDHKRoeIerGo6HI3TKc2wXT4F5jyU0S0SGpyBaQmlwNmD77sSnItZfLUOLHc0uyL5omRn4maQzobCfIVISQ3GkRwRIhhGKPBYojxDCazEB9S/TYwH/gRpqw9MB37A1LGDkuUyDrxXgnKCoGcQzBl7oEpYwdMUrlYLNWLg0RjNBp0AcwBqoogI/iDtNhDiCvGGIZRAhZDjGcJikRFh1ssIuXgElHBZd7yofARkZ8IfjUPB6yVinKAqtayD8KUfUB0vTbBMuzPNlA1qR9Qv43LXiW9RoPENiyEBNw/iGEYubAYMgD0hZl95qy2UmVVp9vHd4YU3VJ4iMxpWyydqrP2WLo8U1UXiSI6+QcD/nT+/H/z+ehBUZYQPiSAkHPYlgKzItWLFyXydDvU88id+6gkRhBCnBZjGMbosBjyMdw1UbssiAgSPW2vR3l8F5j3fy8m3ZOnqDYkmESKi4am2i2nfUW3BEgA0f+gSLgLC6HqsBBiGMYXYDHEeJ+Y1qjoOVE0Z8TZdEuU51wxQP+t58+VwIQKS/qLpseb/IDIppboD4mfsATFyvWVFkEER4Ts4dJ5hmG0DIshxrvRISs0/qJpf3H2gtunSpPEirLz4qgECIoA/OQZc/Uqggguna8ZLptnGMbQYqisrAxbt25Bbk4OmjRtitat2yiyTVpaGvbv3wc/Pz+0atUacXFxHroHvokigqg2qCM0iR86eWDoq9aEkFb8QXqsGGOjNMMwuhZDJGYmT35WnG/arDm+/HIOeve+GA+Of9jlbSRJwuuvzcKRo0fQsmUrlJSU4J23Z2Ps2HEYee0oL90z32i+6HFBpAMBZIWFkD2cFmMYRk+oKobmzv0f/Pz88cqMWQgKCsKRI4fx5BOPo2ev3ujZs5dL25AY6tOnLx7/vydhOj9vauU/f+Ott95E74svRsOGCdAbFAHIzi6CFtGDINKqCNJSREhrRml3+weFB9feVyqilm0ZhvE9VBoYBVRUVGDt2jW4sn9/IWqI5s2T0aZNW6xZvcrlbcxmMy65tI9NCBGdO3cRIik1NRVGxZUvVT2IDXePS6tCiJ4vV/xBLIQsIoiFEMMwSqLaz6PMzEwUFRUhKdG+J0xiUhIOHTyo2DbExo0bhEhq1qx5rT4kOlkpLCyEkVF6TpmWIkSeFmdGiQbpMSLE/iDG1wkLCqgz8lkrpW5sa2BUE0PFRZa0T2io/RdDWGgYiooKFdvm+PHjmDPnc+EXql/f8Zfkd4u+xYIF39guVxZGjHMiRA1R5K3olK8JIb1MnOe0GMMwuhRDgUGWDz6K9FSmsKgQgedTYO5uc+rUKUx/YQp69OyFceNuqfV4Ro+5DsNHjLRdzsvLQ0KC/vxFWsBbUSJvp+eMIoQ4GsQwDKMRMRQbGwd/f3+kpafZLU9LO42EhEZub3P69ClMfv5ZdOjQEQ89NEGkyWojICBAnKycO3cORkfpVJkjoaKkMFLDn6R1EUSwEKoZNkozDKNpAzWJmou6dceqlf8IczORlZmJXTt32lWS7dixHStX/uPUNqdPnz4vhDrg4QmPiF5DvoCaJmo5RmZnDM1Vt/G0GVqvQohEkBaFEKXFPOkPorRYbakxEkEshBiGkYuq9aW33nobJj3zNF5+aTpatW6Nf/7+G23atsVll11uW4fK4vfv349+/S6TtQ31FZoy+TkR2WnZqjV+++1X2746deqExo2bQIvQF1pGTrGhokN6qz5TWwT5ij+ITdIMw2gNVcVQYmISZs9+G3///RdycnKEb4dETeVITqdOndGgQUPZ21DEqEfPnuL8qVMn7W6vRYsW0Cta7jVkJNwVQQSnxRzDQohhGC1ikqz5JsYOMlBHRkbi23WHEBp24Qsyvah69OZknv1EdSIlt/oysX2uY0FTV2RIrhjKPnMWzuLt6JDWUEsEGaVs3t1qMS6ZZ/RORLDZK99Ju49mIDwiwuX95OfloX2zOOTm5iLCjf0YDW7Dyvg0SoggrUWDCE6LMQzDyIfFkIegX9uOokNaRA3vkNroJRpEaNEkLW7Pwybp2mCDNMMwSsFiyMvQl1VtqTI18RVBpGY0yFfSYuwNYhhGT7AYMqCJmr6kXfENGV0Q6UkEEZwWqw5HgxiG8QQshhjDCyKlRJAWhZCWokEEp8UYhtEjLIachL4Qaqoo01qvIXeiQ0YQREoKIC2KIF8ySXM0iGEYT8NiiDGUINKKCBLbcjSIRRDDMLqAxZAC0C/rmnoNGaH5olVcaFkUKS2AtCqCCI4GMQzDKA+LIQPjbqpMy1EiTwggvafEtOIN4pQYwzB6g8WQxsrrlZ5RprQgItQSRZ4UQARHg7hvEMMwvgmLIR3ibKpMSUFUkyjxlDjytPgxgggiOBrEMAzjHiyGPIjeulArIVpcFUbeEj5qiSBf9QZxF2mG8S3Wr1+PDz74ADt27EBsbCxGjRqFe+65x24A+6lTpzBr1ixs2LABpaWl6NChA5544gnxXy1YDPkISkeHtCRqvCmCjBQNclcIsTeIYZjKrFy5Ek8++STuv/9+jB8/Hvv378cjjzwihNG7774r1qHZ8AMGDED9+vXx8ssvIyQkBG+++Sb69u2LPXv2oGHDhlADFkN1EBMcgBJ4Fzm+IVeqyrwliLSKt0WQWtEgtcdpcDSIYXyT3r17Y+3atbbLPXr0QF5enhBEJHgCAgJw/PhxIXpWrFiByy+/XKz33nvvCXFEkaIRI0aocuwshnQ8o4wFke+IICPMFOPmiQxjbAIDq39+mM1mu8uNGjVCixYt8Ouvv+KKK66AyWTCTz/9hLCwMHTv3h1qwWJIo72GlK4q87UIkbsCyGgiSNweR4MYhjkPRWwqExQUJE5Kkp+fj1dffRXXXnutiAoR9P/PP//E8OHD8emnn4rbpNQZRYoSExNVe35YDHl4JIenTdSuNmE0qiBiEVQdFkEMYxzCQ/zdi7KWWbZt3Lix3eIpU6Zg6tSp1VYnf09qaqrD3ZHH56+//qp+M2VluPHGG4XQoTSYlfLyctx6660IDw/HwoULhWfonXfewbhx47BmzRrExcVBDVgMaThVJjc65OuCSAkB5I1IkC+lxAhOizGMdklJSUFERITtcpCDqNBHH30khI0jrBGfytD6N9xwA/bt2yeEEvmBrFB67O+//8axY8fQpEkTm7eI0mcff/wxJk2aBDVgMeTjWIWE3kSRUgLIVRFk5JQYiyCGMT4khCqLIUeQv8cZzp07h7Fjx2Lbtm1CCFWNQBUUFIj/MTExtmX+/v7iWKzXqQGLIYPg7swyPYgiJQWQ2B+LoGpwlRjDMK5iFUJbtmwRQsga+anMpZdeKlJjVFb/4osvCoP1119/jUOHDqF///5QCxZDKpuolUqVKTXEVUuiSGnxo0URZISUGKfDGIYhli5dikWLFqFBgwYYNGgQKvP7778LcUSRom+++QYTJkwQzRmpAo0aL86ePRsDBw6EWrAYkgF9kWQVlfpUJ+rKQsSbwshTAshbIsiXUmIsghiGqWq2ph5CNZGQkGA7T72E6JSZmSmEEF1HJfZqwmLICxVl7uLt6JC3okWeFD7uCiCjiyCCU2IMwyhJeHg42rZtK3t9GtehFVgM6aQBo9qCyFviRSlYBDmGRRDDMIw9LIa85BvydqrMKgY8IYq0ijcFEMGRIIZhGGPAYsig0SFPR4mMIIA8mQojOB3GMAyjD1gMuWGidtY3pER0yFVBRBhFFLkrgAgWQdwwkWEYxgqLIQ3NKZM7vNXVuWV6FUVKiB9Pp8KMEAni6jCGYXwVFkNeRinvkDuDXLUuipQSP1ZYBPHoDIZhmNpgMaQx5EaHlJhsryVRpAUBpNVIEMG9ghiGYTwHiyEP+YZqS5XVFR3ypiByJEQ8JZCUFj1aFUFKdI0mWAQxDMN4HhZDGsVZQUS4K4q8KVp8VQCJ22URxDAMoylYDCmAJ6JDzgoipaJERhY/zgogb4sgT3aMJtggzTAMUzMshjw8p0wNQUQYSRS5I360LoA8mQqrDAshhmEYx7AYUnlWmVxBRLgiivQmjNwVPp4WQFqIAhEcCWIYhlEOFkN1UNOXkqPokCvpMmfK7Z2NEjkSGFoRR0oJH1fEjxoCiOBUGMMwjPZgMeQl1BZEaoojJUWPuwJI6TSYFqJAnAJjGIZxDxZDLuJsdEgOzggiwl1R5GmxokcBpJcoEAsghmEY5WAx5ABJksT/goJ8hIUBBUXnqq1ztrisxm2LHIihaDNwKr92sRPnD6Tm1bzfqkSd/77MyFO/aaKniIuwFz+lhQWyt02MuCAois+W1LpuQvgF4VJU4Hjd2PMi6GyB4+co5ryQKcgvrlUAFeRX3zYsyLJtfmnNz2n4+e1lvkQYhvEWpWaEh4fDZDJ59Gby8/JU3d6osBhyQHZ2tvjfv3sbbz4fDMMwjE7JzMxE/fr1PbLvwMBANGzYEO1bNXV7X7Qf2h9zARZDDoiNjcXVgwfis8//h9DQUOiRwsJC3HfvXfjo489Qr1496BUj3A8j3AeC74d24OdCe89FcLDnLAe07yNHjqC01P3ZliSEPHmseoTFkAPMZrP40oqMjNTtl5e/vz8CAgIQERGh2/tglPthhPtA8P3QDvxcaO+58HSKjAQMixjPYPbQfhmGYRiGYXQBiyGGYRiGYXwaFkMOoJDnDTeMFf/1ihHug1HuhxHuA8H3Qzvwc6EdjPJc+DImyVpDzjAMwzAM44NwZIhhGIZhGJ+GxRDDMAzDMD4NiyGGYRiGYXwan+0zVFJSgpMnUxERESmrY+jBgwdQVmY/A4G2i49vADXJz89HRno6YurXR1RUlOwGYXTfIyOjEBcXB7Uh29qpU6dw7tw5NGjQAEFBQbWuf+LECeTn27eUp949TZs2gxY4fvw4zp4tQKtWrUX/kdqg11RKSopogpaUlAStQMdFr3kaL5CU1NjhelmZmUjPSK/Wo6tNm7ZQC1ffq/Q6pOeiorwcjZs0gZ+fH9Ti7NmzOH78WLXlycktHL4/XNnGW9B7m9631MBW7mdOWlqaGIeUmJikmd46eXl5yMzIQKPExFqPqby8HPv376u2vFGjRNG7jtEePieGcnJyMG/u/7Bu3To0aBCPjIxMJCY2woSJj4kW5Y547dVZ8PP3Q1TkBcHRt18/DBkyDGpw+vQpfP7ZpzhwYD/q148V4qZN27Z49NHHER4e4XC7X3/5GXPmfI64uHhkZmagc5euYhu1Piz/+fsvfP3NVwgMCEBFRYV4fq6/YSxGjBjpcJuvv5qLvXv32j1fzZsn455774PaHDp0EM89+4zoEvvJp18gOjra4brbt2/Dm2+8juCQYBQVForn8ZlJzyI2Vn2BSq+tZct+Q/cePfH005Mcrvf3P39h0bcL0axZc9syei1NnjINauHKe5XePzNeeVl8+VJFUEWFhP974knVRB0JuhemTUHbtu3slj/y6OMOxYQr23iDFSv+wJwvPkdEBM3tMiMxMRETJj6KkJAQhz/WXp01Q4iJmJgYnDlzBvfe9wCuuOJKqEVRURHef/89bNq4Qfxoyc3NFZ9TAwdeVeP6xcVF4nOA3heVRdP1N9yIrl0v8uKRM7KRfIwDB/ZLf/yxXDp37py4XFRUJE2Z/Jz0zNNP1Lrd/ffdIy37/TdJK2zbtlXavn2b7XJeXp704IP3Sx+8/16t933M6JHSurVrxOUzZ85I995zpzRnzueSWvz001JxHFZWrVopjR41Qjp06JDDbWbNfKXW+6kWhYVnpfEPPiAeT7oPle9XVQoK8qXbbh0nzZ37P3G5tLRUevbZp6XJzz8rqc3atWukRx+ZIE1/Yar0yisv1bruokULpSf+71FJS7jyXn38sUekl1+abvtcoNfXPXffIZWUlEhqsHXrFvEa8vQ2nmbdurXSdWOutX3mEBs3bpDS0tIcbvPf996RHnroP1J+fr64TM/l9deNklJTT0hq8dJL06VHHnlYys7OFpfpdbF8+bJa39/0XNBnLqMPfM4z1LJlK/TvP8AWAifVfvkVV+LAgQMitFlXSop+fVH0Qm06d+6CTp062y5TOoPSRFlnzjjc5s8Vf6Bx4yboffEl4jJFLQYOHCSWq9VhYejQYXbRk44dO4n/2bXcD2uak56L9PR01Y69Kh9+8D569uqFTh0vPC+OWL9uHYqLizFq1BhxmaIRo64dg507dyA9PQ1qkZGRgU8+/hCPPPIY/GX2TCkvr8DRI0dEGoTSIVrAmffq4cOHceTIYYy57nrb58J1190ghjVv2fIv1CQ19QSOHj0qXu+e3MZTfPPNV+jTp6/tM4fo0aMn4uPja1yf0psrV/4jPhfCwsLEsgEDrxKppb/++hNqQNH3zZs24q677rFZESitPWDAwDq3pTQyRYspbc5oG59Lk9UEfWjGxsXV6RFYtGghVq1aKVJUzZOT8dBDE2tNrXmDffv2ipTMgf37sWf3Ljz51DMO16UP/OQWLeyWtWzVCnnz85CVlalaeoZCzpSmKCgoEGm8Dh07onOXLrVus3r1KvGBn5WVhbDwMDz44MPo0KED1GLFH8tx4kQKxj80ATt37KhzfXouEhIa2c0po+dCXHf4sCpeNPox8OYbr2HU6DFo0lT+ZGzyqbw5+3UUni1EaWkJ7rjzblx5ZX+oiTPvVXourN4aK/VjY4VIp+t6974YavHiiy/Az+wnUtrDh4/EuJtvqXP+lSvbeOp9ffzYMYwZc514b6elnUZs/VhE1uJtTE1NFSKuRXJL2zI69hYtWor3hRrs2L5dpPQ6dOgovI30Gm/YMEGWteCDD94XAoo+3y6+5FLcd98Duh38bXR8XgzRL/Flv/+GBx54sNYH6oYbbsRll18hDLFk3p058xW8/toszJj5qqpGy/nffC1+BdOvwX6XXW7n3ahKfkEBWoeF2y2jiBJRkF+gmhiiL5yFC+YjJzcHxUVFuOee+2rt5Hppn754cPxDCA0NE78kKZIxa+YreOvtd2WbyJWEHvsvv5yDF6a/JLsDLX05WB97K9ZfwvQ8qfUrnj70hw69RvY2JCDefe8DYXwnfv75J/z3vXeQ0DABbdvZe1e8hbPv1YL8fCFKq5rdw8LDxftCDeh1/NLLM2z+H/qcenH6NCHQhg67RrFtPIk1urtn927hQYuOjhGigDwz5BmqyYBMni3rY18Zeq9Q5FENyLNEkalZs2YIcUevk+zsM7j99jsx8KpBNW7j5+cvoqv0mUyQiJo65Xl89unHeHjCI16+B4wcfC5NVhkKX86c8QqGDRuO/nWEPOl664clGZRvvvlWHD58CKdOnYSakFH11dfewAcffoKDBw7gvXffdriuv58/yspK7ZZRVElcF6CeLqYPR/oQf++9D3D/Aw/ijTdew7ZtWx2uT2F3EkIEiY8777oHhYVnsXXrFqgBPeYU+ieBs2fPbqSkHBfLDx7Y7zDlRV/KVZ8LawVUXRVonhKkS378QRhC9+7dI+4HfTHRic5bXyc1PXdWIURQeoOqf1avWQW1cPa96ufvX+P9o2VqPBcEpbwrG6Epfdynbz8R7VJyG09CJnZiz549ePud9/Da62/inXffF8bo+fO/dvgZRVR9b5So+FzQ/Th9+jRaJLfAf9//UNyX226/Ex9++L7tvV4VEnpWIUQkJCRg5LWjsGbNalEowmgPnxVDhw8dEpUX/fv3x+133On09lSWTlCaRgtERETgyv4DavU4UEVJVU/RmawsEYamSiYt0KtXb5E+2rpFvrChD57g4BCcOaPOcxEVFY2TJ09i7pf/EyeqnrGmajZv3lTjNlTNl5VV5bk4f/xqVP6I1ESLlli6dIntfqQcTxEnOk/RR7lQhIJeV1qhrvcqPd7kdaKy6copw9ycHJE+1wpUHUdRCk9voxQUaabPlr59+9l+vFCLgx49eoloUU1YX/tnanhvqPVcWP1NgwZfbVtGvlNqIbFv716n3hcksOlHE6M9fDJNRr+Cp02bIkLpFFWoCYqyRERGCO8GGV2rhnS3bdsi3gy19WDxJDUd0+lTpxBWKQ1Gpj3qeWPtM0Jl9N98Pc9u2w0b1qN169YOy1w9CUVC6MOy8i8+OracnGy7FBJ5cejLiX750pcWGaYrp6PolyZFhpo0ke9zUZKqPq0t//6LF1+chqeeftZmDrf2HbH2GenSpSsWLPhGiHKrj4ueC0rXUH8ib0MRBYrOVWbGjJfF/8ql9eRDIVFhLTmv+jqktBRFYSjaqgZy36uV39/t23cQr8GNGzfYTLHk+6Jy6i51eNe8dT8omrBjxzY0adLE4ftbzjbehD5T6HWVVeVHCvkTIyr12qn8/iavFkUW6bm4qFs3cT2J0v379mHggJrL2D1Nly4X2QSZtUdQXm6u+Cyy3o+q7+8aX4dbt4gfTlXT44w28DkxRDnraVMniw+ISy/tI1IAlSvNrF+y5DOgXzQUNSKTMv3Kv/zyK0XfCzIrL168CMNHjJTVsNETfPTh++IXL/ky6IOcPrx/++0X/Oc/423r7N+3X3wpv/X2e6I3BqVAyKA8Y8ZLwhdC92P9+nV4fvJUVe4DeTVeeukF9B9wFRIbJSIvLxe//PKT+BClChIrc+d+icKzZ4Un52xBAaZOm4yrBg4Sjc8o9fHtwoXo1LkzunXrDq1CDfGo7wh5JS6//ArxvFFq7c03X8PYm25Gfl6e8H+R2ZUqVbTKX3/+KV77876aLy7Te6l79x7CpEy/eH/4frEQ5EOGqtN/S+57tfL7m76crh01GnO++AxSRYUQFuQBo742VH2pBuSDI6Nt+w4dqf8J/vhjuTAXUxrZ0ftbzjbe5uZbbhW+JXrsSbTt2LFdpLOnTpte4/vbus1rr84Uzx81v/zh++/E89C332Wq3Afx2XnVILz91mzhR6O06uLvFgl/5kUXdavx/f3777+JwpyePXqiXmgoNm/ahD//XIEHxz+sipmdqRufm1pPb0Qy6zr6hW9V/mTIpS9Ya6M2Kq/8Y/lypGekiYqIPn36oUvXrlAL+lVCFUzkraHQa4MGDTFg4EDRfNAK/VKhZmeVm67Rr6zvFi/CsaNHxa+aq68eIn4ZqwV1mSWBduz4UYTWCxWC9KpBg+2qrObN/VL8Src2VSQfzi8//yyqmCg92KlzF/HFRb/+tYD1cX/6mWdtvwLpw/Lll6bbNV2j523pkh+F0ZUEEPk7+qn0gV8T9Lhbv5ys/PXnCvzzz9+2pooUnfjll5/FL3e6D/T8XT1kqKodg+W8V6u+v+ljkFpMrFu3VvzK73rRReI6tXwqld/fdJ5EwbCh1yA6Jsbh+1vONmo9H7/8/JNoVRDfoAEGDx6C5ORkh+9vYtvWrVi2/HfhWaP0LbWgsBYYqAG9JpYt+12U2JvMZrRu1RrDrrkGISH1HL6/KdK7du0a8UOnQcOGuOqqwWjWTBtd8pnq+JwYYhiGYRiGqYw2fkozDMMwDMOoBIshhmEYhmF8GhZDDMMwDMP4NCyGGIZhGIbxaVgMMQzDMAzj07AYYhiGYRjGp2ExxDAMwzCMT8NiiGE0xt49e8SsM+LYsaOia6/aUBM/ag65etVK0SjTKKxft1Y0zGMYxrdhMcQwGoJmf82c+TKCg4PE5dWrVznsmO6Io0ePCuGiZPfdSc88ic8+/USMb6FjNAo7duwQ8/oYhvFtfG42GcNomQXzvxEz82JiLHO0aHilsyMhVq78G4cOHUTHjp0UOSba1+HDhzF33jeqjtnwBNeOGoXxDz6AESNH2UbWMAzje3BkiGE8AM1hWrXqH5SVldmW0ayi9PR0h9vQoNOVK//BlVcOsBsS2a5de9tlEiW7d+8Sc81oHtiWf/8Vk+KtpKaeQOqJE8jNyRW3T6eszEy7CeEU3aF5UZT6qox13yUlJWK+Fa13/NgxrF+3Dn5+fti0aYPYH01Dd3V/laNNNCR544b1YnhyVWhq+enTp5GbmyvmCdJ+qu7fuh8azkrr01Txmq6v7XZiY+PE47ts2W81PCMMw/gKHBliGA9gNpkw+803MOnZ59GtW3ebcFj2+2949rnJNW6zY/s2EQVKbtHCtozSZOQh6tSps21Q6rbtW2GCSUz1JgGVlnZaDE6lgZanTp3CqVMnxfIN69eLbeLjGyA8IgJvzX4De/fuRYuWLZGeliamZz8z6VlxvW3f27bQwHMxZZyGanbt0hX79u8VAsi6v4svvlQIE1f217v3xeJxeOWVl8Rg1/i4eCGkaLjlQw9PtA3b/fSTj5GQkCD8SUmNG+PokcOIiorGSy/PQEBAgFgnJeU4Zs54xTKUtHFjnEg9IQZ6Dho02PZ413U7RMdOnbB2zRqMG3eLB14JDMPoARZDDOMBIqOihDihiIRVDA0cOAhfzZsrIhSNGiVW24YiKYmJiXZf1DVx6uRJvPzKTDEh3jqBfdGib/Hkk0+jR4+eIhJCqa3HHn/Cts2X/5sjoiz/ff9DBAUFiSnt/33vHXzy8UdCsFlJTU3FC9NfQvv2HWzLAgID8dqrMxXZH4mqV2fNwJVX9Md1199gm3z/f48/KiauD7xqkG3bnNwcvP7GbISEhIh1HvzP/UIcXnHFlSLiQ0KoWfPmeOSRx4SIJFG0ffs2p2+naZOm+Obrr0QUzyq0GIbxLVgMMYyH6N69h0i/3CdJImoSHh4uvrRPnDhRoxjKy89DaGhYnftt2bKlTQgR7Tt0xPJlv9e6zYoVf+CSSy7Fv5s3QQKEeKFIy5o1q8V5Oj6iWbPmdkJI6f3t379P3P+Y+vWxltY9v23DhgnC9F1ZpFx+2RVCCBH0uDRt1syW6qLUGEXAnnt+ss1TRf+twtOZ2wkNCxPXFeTnIzomps77zjCM8WAxxDAeolv37pg//2scOngQLVu1wr69e0X0oUEDSxqpKmROJn9NXYSFhdtdDgjwR1lZqcP1aZ95ebkibVRQkF9NsFFExRoRiY6OrvP23dlfRnq6iHxt3fKv3XISipQOs7+f9sIwwD8AZaUWD1ZmRobwMTVo0LDGY3TmdkqKLY958HnhxTCM78FiiGE8BKXJKFpC5l7yAc2d+z+0atVKVIjVRKNGjbBm9WrFj4OECUVN+va7zOancYQ1ouOp/YXUqydSWHfdfS+ioqLgKhQpolRZYeHZGqNpztxOWnqa8F9Zo1AMw/geXE3GMB6ChEC3bt2wYeN6LFjwDQ4ePIB773vA4fqdOnURFVG1VZzJISQ4xBZBIShC0qVLV5FKIwFRmays6hVYdeHO/tq1ayciYL///qvdchIuVIEnl7bt2gqv0l9//Wm3nHxMzt7Ovr170LlzF9m3zTCM8eDIEMN4kG7dewh/DZWoPzj+YREtcgRFhjp06CjK10ePvs7l26TbWLRoIX76aQkiIyPRrm173HnXPZj8/CQ89+wz6NfvMlRIEnbv2ikEw4SJjzp9G67uj6I4JAjf/++7ogKtbdt2QpxQ2f2o0WPQp09fWbdP+6Goz8cffSC6dTdv1hwHDh6An9kP993/gOzboZTfpk0b8cyk55x+DBiGMQ4shhjGg3Tt2hV9+vZDmzZtMGDAwDrXv/6GG/Heu+9g+PCRIh1VtelicnKySOlUhkzB3bv3tF2+qFs3PPjgQ9i1axf2Fe4Vpe6tW7fBm7PfwZ9//oHDhw8hNDQUV17ZHz179a513wQ1gKRGkJWhsndX90fVYMnNk0VPJap8i4uPFwKqSZMmtnWoKq5BQ3s/UIeOHRFbP9Z2eeDAq9CsWTOsWrlSGKZbt2krjsGZ2/nrrxVo3jxZlmmcYRjjYpKojIJhGM3wzTdf4aKLuqFNm7ZqH4rh+eLzz3Bl//4OfVwMw/gGLIYYhmEYhvFp2EDNMAzDMIxPw2KIYRiGYRifhsUQwzAMwzA+DYshhmEYhmF8GhZDDMMwDMP4NCyGGIZhGIbxaVgMMQzDMAzj07AYYhiGYRjGp2ExxDAMwzCMT8NiiGEYhmEYn4bFEMMwDMMw8GX+H6LEOyFa0b/cAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nus = np.linspace(2.5, 6.5, 41)\n", "ps = np.linspace(0.005, 0.12, 41)\n", "surface = np.array([[loglik_fast(mz, None, nu=a, p=b, kernel=kernel) for a in nus] for b in ps])\n", "rel = surface - ll_hat\n", "\n", "ink, ink2, grid, accent = \"#0b0b0b\", \"#52514e\", \"#e1e0d9\", \"#2a78d6\"\n", "fig, ax = plt.subplots(figsize=(6.4, 4.6))\n", "levels = np.arange(-30, 0.1, 2.0)\n", "filled = ax.contourf(nus, ps, np.maximum(rel, -30), levels=levels, cmap=\"Blues\")\n", "ax.contour(nus, ps, rel, levels=[-chi2.ppf(0.95, 2) / 2], colors=ink, linewidths=1.5, linestyles=\"solid\")\n", "ax.plot(nu_hat, p_hat, \"o\", ms=8, color=accent, mec=\"white\", mew=2, label=\"MLE\")\n", "ax.plot(nu_true, p_true, \"x\", ms=9, mew=2, color=ink, label=\"truth\")\n", "cbar = fig.colorbar(filled, ax=ax)\n", "cbar.set_label(\"log-likelihood minus maximum\", color=ink2)\n", "ax.set_xlabel(\"ν (interference)\", color=ink2)\n", "ax.set_ylabel(\"p (escape fraction)\", color=ink2)\n", "ax.set_title(\"Log-likelihood surface, 2,000 simulated meioses\", color=ink, loc=\"left\")\n", "ax.legend(frameon=False, loc=\"upper right\")\n", "for side in (\"top\", \"right\"):\n", " ax.spines[side].set_visible(False)\n", "ax.tick_params(colors=ink2)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "394d0e41", "metadata": {}, "source": [ "## Profile likelihood for $\\nu$\n", "\n", "For each $\\nu$, maximise over $p$. The interval where the profile is within\n", "$\\chi^2_{1, 0.95}/2 \\approx 1.92$ of its maximum is a 95% profile-likelihood interval. It does\n", "not rely on the quadratic approximation behind the Wald interval, so the two can be compared." ] }, { "cell_type": "code", "execution_count": 7, "id": "f012472e", "metadata": { "execution": { "iopub.execute_input": "2026-10-26T15:47:59.259311Z", "iopub.status.busy": "2026-10-26T15:47:59.259250Z", "iopub.status.idle": "2026-10-26T15:48:00.364387Z", "shell.execute_reply": "2026-10-26T15:48:00.363970Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "95% profile interval for nu: about [3.40, 4.30] (grid step 0.1)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkoAAAFqCAYAAADlZ3UyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAb6FJREFUeJzt3QWUlGUXB/D/dnfvssvS3SldgkqIgKCi2IiJAordKCiCoGIriEh8NCithHRId7Owy3Z3fOc+wwwzsJ2zM//fOXP2nX733dmZO89zn3st8vLy8kBEREREt7G8/SIiIiIiYqBEREREVAiOKBEREREVgIESERERUQEYKBEREREVgIESERERUQEYKBEREREVgIESERERUQEYKFVjOTk5SE9Ph9QMzc7ORkZGhsH1+V1WFsV9jvJ8Xvnd5HeU37U4l+Xm5sKYyD7J8SgJOXYlvQ8REVUMBkrlJDMzU30oykk+6Cqy4PnWrVvRvn17+Pj4ICgoCCdPnsSbb76JwMBAg9t9+OGH8PX1Lbfnze/x8rts3LhxqF27drk85+XLl+Ht7Y25c+cWetmRI0fUZUuWLIGxkGBH9mnSpElF3jY2NhajRo1Sf0/5O7711luVso9ERFQ46yKup2Lq3r07jh49CltbW3U+KysLoaGhGDFiBCZMmAA7O7tyOZapqal49NFHMWjQIOzYsQPW1po/oY2NDezt7Sv97yXPX9nPa2FhoY6n9nc3BRJM/fvvv9izZ48KloiIyDhwRKkc1apVCzExMeoUFhamAppPP/0UL7zwQrk9x4kTJxAVFYUBAwYYBAryQXv16lVUtnfffRfXr1+v1OcMCQlRx/jhhx+GqZCgV0YJGSQRERkXBkoVxNXVFePHj8cdd9yBhQsXIi4uLt88Ghl5ym+aLr/LZHrv2rVrmj+cpaVuqk9/yq+k9PN8SqMk+Uhy21tzi7RKkluUX45Sfor7mMU9BsW5XUmnXGUf5XeRwE9GBbV/y/wep6DnL+h1VdDvn1/elLy25D5ERGSIgVIFq1u3rvrAioiIwM6dO1XOyrp16/D111+jQYMG8PDwwH///aduu379evTq1UtdJqeOHTti0aJFuscaNmyYymPRbsvog/bk5+d3W45SQQ4dOoT7778f/v7+6nmaNGmC6dOnlyqvqrh5UHv37kX9+vVx7733IiEhQTeN+N5776Fx48Zwd3dXv4OMEkkOUmHyy1HSN2fOHDRq1Ahubm5o06YNNm3adNttZBRszJgxanRKnltGA19++WXEx8eX6nayT8OHD4eXl5c6Hk899RSSkpKKPC7Lly9Xfz8JgFesWKH7e0reWXGfXz8/68cff1TH09PTU03l3UqCIXlNal9HWp988gnq1KmjRiuJiOgmBkoVTD7wrKysDIKYX375BSkpKdi1axe2bNkCR0dH9SE5dOhQNG3aVN3nwoUL6N+/P5544gl8++236n4rV67EH3/8obbl9tppPjk988wzxdofec7evXurvCLZjo6OxtSpUzFt2jS88cYbFXIM5AP87rvvRp8+fbBq1Sr1IS4jGBI0LViwAF9++SUiIyNVfo4EF3feeaf6nUpj6dKluHLlCrZv366OoeSJSfAlydJa8hz9+vVT012LFy9Wx2D27Nn4888/1ZSmdoSsuLeTwO+uu+5Szyd/T3n+++67T40oFmXIkCHqd5XXh9xH+/eUQK+4z68lI5cy/bpt2zYVlEtgdSsZtZLHlOBRRpa0XnvtNTViJQE8ERHpyaNy0bFjx7ymTZvmpaWlqVNYWFjee++9l+fk5JT3wgsvqNts375dnR88ePBt92/WrFle27Zt83Jzcw0uv/fee/MCAgLyUlJS1Pm//vpLPcaWLVsMbvfqq6/meXp6Glz2zjvv5Lm6uhpc1qlTp7xWrVrlZWZmGlz+9ddfq9teu3atwN8xv8fL77IXX3wxLzAwUG1PmTIlz8XFJW/y5MkGt/n555/V7/H3338bXB4dHZ3n5+eXN2nSJHX+4sWL6na//vqr7jb5XXbo0CF12aBBgwwe7+TJk+ryH3/8UXfZl19+qS7btm2bwW1XrlypLp89e3aJbjd9+nR1/uDBg7cdU7n8/fffzytKvXr18kaNGmVwWXGfX/u79+nTJ684li1bpm6/Zs0ag8tHjBiR16JFi2I9BhGRueCIUjm6ePGibupEEnM3btyopjRkWkufjD7okxGI8+fPqxEkWdGlT1a3JSYm4uDBg2XeP5n+k2k3eX55HpmGkZOM7kgulYwoyChTeZDpxtGjR+Ozzz7Dr7/+iokTJxpcL9OPksfVpUsXlS+j3Q9nZ2c1miKjKKUho1H6ZLpPkt4vXbqku0xGfWTKUZ5b3z333KNuK+UXSnI7+RkQEIAWLVrcdruyKO7za8moXXGPkYworl692uByJycnTr0REd3CdNZXGwHJHzl8+HCRt5MPVX3afBOpi3Qr7WWSDF5W2vyTWbNm4fvvv7/tellyL0FZeZBpobVr16q8lx49euS7LzK1JHlS+WnevHmpnvfWfCkJCCUA0OZFaY93fsdapkhlWlB7rIt7O/lZ2N+utIr7/AW9rgoix0Ny4dasWaPy0uQYaYPkZs2alWmfiYhMDUeUqoDkieiTBGAheTr5jQLp36YsJOFXSF0n/fwm/ZOUNCgPDg4O6oNYHlNyYsLDw2/bFwkCCtqPf/75p1TPe+uIXH7kWOZ3rGVkS55be6zLeruylk0o7vNrlaSu1MCBA9X+SZK9mDdvnhp1e+mll8q0z0REpoaBkhGQRF6ZIpJk7VuXdC9btkyNHrRs2bLMzyMjDq1bt1aJ4DLNVdFkNd2GDRvU6JJM98j0opZMM0oQsHnzZlQ2GeGS0ZpbgzE5LjKyIqMtJbldz549VUC7b98+g9vJ37My9rM0ZPpORqYkMVxqfklld1mpV9bpQiIiU8NAyUhIYUoJJJ5++mmcOXNGfbt//fXXVR7KRx99VG7Vr2fOnKlWRslSdplqkdVgZ8+eVauqZFWa/hRVeU1HSrAkK/skWJLq5WLkyJHq/OOPP66W+UuelgROkpsky99lerCiyHNKHpQsu5fVXzINKPk60npFSjJI6YWS3O6xxx5DvXr18OSTT6ol+fJ7yGq+06dPV8p+loaMRnXq1EkFXfJakABaVj8SEZEhBkrlRPJ7tO1LCiJFIuV28k3+VjI9Jd/uZQm4LN/v3LkzDhw4oOoo6U+HyX3lMeSx9OXXwiS/9iIyMiXBSHBwsArKZNRHairJNJkEa1J7qCD5PV5+l8m+6LdskTwkSd6W/m+DBw9WdaPkfhKcybJ0KZcgyeRykrpMkicjQUJB7Uryu6ywYyuX6093StAm+yOJ8mPHjlXP9/bbb6vjLIGD9rbFvZ3k/Mjxa9u2rQoA5feQ5fmSyH7rcxdEjuGtr5/iPn9hv3thpDyDBKiSVD9//nyTaglDRFReLGTpW7k9GhEREZEJ4YgSERERUQEYKBEREREVgIESERERUQEYKBEREREVgIESERERUQEYKBEREREVgIFSGUl1BemPxioLREREpoeBUhlJY1cp0ig/iSh/mRnpuHzptPpJRFSdMFAiIiIiKgADJSIiIqICMFAiIiIiKgADJSIiIqICMFAiIiIiKgADJSIiIqICMFAiIiIiKgADJSIiIqICMFAiIiIiKgADJSIiIqICMFAiIiIiKgADJSIiIqICMFAiIiIiKgADJQDZ2dlISUkp6BgRERGRmbKGGcvKysJPP/2ALZv/Uee9vX3w7HMvoEmTJlW9a0RERGQEzHpEaf78eTiwfz+mTZ+B3+ctQIeOHfHpJx8hPj6+qneNiIiIjIClOU+3bVi/DgMHDkJgYBCsra3xwAMPwcLCEpv/+buqd4+IiIiMgNlOvV29GobU1FQ0bNRId5mNjQ3q1a+HM2dOV+m+EZVVamYu4lNyEJeajdjkHMSn5iDuxvms7Dw0DrJHm1BH+LnZ8GATERXCbAOlxMRE9dPV1dXgcjdXN0RFRxWa1yQnLQm2iCrThagM7D2fipjk7JsB0I0gSH7KZelZecV6rGBPG7QOdUTbWo4qcArxsoGFhUWF/w5ERNWF2QZKWjk5Obedt7QseEZy6ZLFWLRoge68ftBEVFHy8vKw/UwK5m6Pxb+ny2+F5pXYLFyJTcCKAwnqvLeLlQqY1KmWI+r72cHSkoETEZkvsw2UvLy81M/4+AQEBdXQXR6fEA9PD88C7zdk6DAMHHSvwchUQEBABe8tmau0zFysOpigAqTzkZmF3tbKEnB3tIKnkzXcnazg4WgFD/npZK22tZfl5gEHLqXiwMU0HL6Shszsm6NP0Uk5WHckSZ2Eq70lWt0InGTUqVGgPWytGTgRkfkw20DJ3z8AHh4eOHzooK4cgEyjnT51Ck888VSB95M8JjnpJ4UTlbfrCVmYvysOi3bHIyHNcNQzyMMGwzu4o46v3Y1gyFoFRC72lsWeNuvawFn9zMjKxdGr6dh/IRX7L6biv0tpSMnI1d0uMT0XW04mq5NwsLXA0LbuGN3TG17OZvv2QURmxGzf6WR6bcjQ+/H73DkIqVkTNWrUwPz5f8Dd3QPduveo6t0jM3XkSpoaPVp3JBHZN+MVpU2oA0Z18UTPRi6wKqfpMDsbS91Um8jOycPpiHTsu5CmAic5Sd6TVlpmHn7fEYcl++IxqrMnHu/mBRd7q3LZFyIiY2SRJ8kPZmztmr+wceN6JCenoG7dunhk1GPw8/Mr9v1l6s3NzQ0JCQm3JYYTFYcEJ5uOJ+G37bE4eCnN4DprK+Ce5m54uLMHmgQ5VPoBlbeHC1GZ2CdB04VUbDqWhDS9RHFXB0s81d0bD93hAQfbgnP7MjPSERFxGf7+IbC1s6+kvSciKjuzD5TKioESlfq1k5aDJXvjMW9nLMLjDadwZSpteHt3PNDRA76uxrOEPyopG9//HY3/7Y1Dtt6MoI+LNZ7t7Y0hbd1hY3X7aBcDJSKqrhgolREDJSqpnNw8zNoUjTn/xqipLH11/ezwSGcPDGjpBnsb460HGxabiW82RqtEc/0xaSk38EIfH9zTwtVgtRwDJSKqrhgolREDJSqJlIwcvLrgmi45WqtbAyeM6uKFjnUcq1Udo7PXMzBzfSQ2HTf8fer72+HFO33Qs5Gz+n0YKBFRdcVAqYwYKFFxXYvPwvNzruB0RIZuOf+wdu54pLMnavnYVesDKWUGZqyLxK5zhgVYW4Y4YGw/H7QMsmKOEhFVSwyUyoiBEhXH4ctpeGHuFcQk5+jqE00fWQMd6zqZ1AHcdTYFX66LxJGwdIPLh7RywNMdMpnMTUTVjvEmQRCZiDWHE/Hoj5d0QZK0CfnjuVCTC5KE/E7znwvFzIdroI6vre7yQ1c0q/l+2ByNzFvrHhARGTEGSkQVuLR+1qYoTJh/VVf9ul0tR8x/NrTaT7UVRnKSejdxwbKxtfHJ/QGqQKbW+qNJGPX9JYTHs/UPEVUPnHorI069UX6k4vXbS8Lx1yFN82UxpK0b3rk3wOxagEiQ+OeBKLTxi8cLi3NwLlpT/mDqg0HoWMf0RtWIyLRwRImonEUnZeOJny7rgiRZxDb+bl98OMT8giQhv3P/FppirL6ummYAUu376Z8v4+ctMWrkjYjIWDFQIipHZyLS8eCsizh4WZOT42BjgRkja+CJbl7Vatl/RZkyPAhd6mtGkaQ577S1kXh53lUkpxv2syMiMhYMlIjKydaTyRj53SVVBkD4uVpj7phQla9DGtK4d9ajwXiut7fukGw8loQR31xUNZmIiIwNAyWiMpKpI2lk+/xvV5CSoVnR1STIHgueD0WjQPY1u5U09H2+jw9mPVpDlUkQF6Mz8cCsC2qFIBGRMWGgRFQGWTl5+GhFBCavvq6mkkTfpi6YM7qmUfVoM0bdG7pg0Yu10CBAswJQ2rnICsEpq6+r40pEZAwYKBGVoants7OvYOHueN1lo3t64YsHg+Bgy3+t4gj2tMW8MaEY1MpNd9lv22Px5E+XVANeIqKqxndzolL2bHv0h0vYeTZFnbexssCn9wdibF9fg2awVDQJKqXe0ruD/WFtpbls/8U03P/VBRy4aNgShYiosjFQIipFTtK7S8J1PdukJtDPT4VgUOuboyJUMrIicEQHD8wdHaqS4IWMKD3+4yX8vj2WJQSIqMowUCIqoXk74rD2SJLadrazxNxnaqJNqCOPYzloHuKA/71YC+1ra46ndDv5dPV1vLbwGlIz2fqEiCofAyWiEjh4KRWf/3Vdd/6T+wNNuh1JVfBytsaPT4So2lNaUrxT8pZYb4mIKhsDJaJiik3Oxvj5V9Uoh3i8mydrJFUQaysLVc18+sggONlp3qYOX0nHM7OlBAOLUxJR5WGgRFQMObl5mLjwGiISNCux2oQ64OW+vjx2FaxvU1c1tenuqMnyPngpDWNUsMRpOCKqHAyUiIrh203R2HFjhZuXs5UqASCjHlTxGgTYq2R5NwdNsHTgYhqem3OFOUtEVCkYKBEVYdupZHz3T7TmH8YCquu9D4tJVqqGN4IlVwfNW9a+C6kqWEpjgjcRVTCzDpSys7Oxffu/eP/9d/D8c88gIiKiqneJjIz0bZMpN22D+7H9fNC+tqapK1UuaQfz05MhurYne8+nqrYxDJaIqCKZdaD07bffYOeO7WjTpq0KkiRwItLKzM7FuHlhSEjTJA/3aOSMJ7reXIlFla9JkAN+eCJENdcVu8+l4sW5YUjPYs4SEVUMsw6Uxox5DhNenYhGjRpX9a6QEfr8r0gcCUtX2zU8bFQpAFbdrnrNgh3ww+MhqoaVkOroEixlMFgiogpg1oGSjQ2bllL+/jqUgD92xqltW2sLtUxdm0xMxlGY8vvHg+F4o6fejjMpeOl3BktEVP7MOlAqjaysLKSmphqcyLScvZ6Bd5eG686/OdAPjYMcqnSf6HYtazqqYMnBVrP68N/TKXh5XpiaMiUiKi+apkom4puvZ+L48WOF3mba9Jmwsyt9JeWlSxZj0aIFBoETmQ6pz/PKH2FIy9Rkb9/b2g3D2rlX9W5RAVqHOuL7x0LwzOzL6m+29VQKXvnjKqY/VEONBBIRlZVJBUoPPfQw0jM0jUoLYmtrW6bnGDJ0GAYOuld3PjExEQEBAWV6TDKeZrfvLQ3H+chMdb6+vx3euddfNWwl49WmliO+fTQYz86+grSsPGw+kYzx88PwxYMMloio7EwqUPLw9KyUvCb93CaulDMd83fFYc3hRLUtbTO+HFkDDjdyYMi4tavthG8eDVa1ldKz8vD38WS8uuCqqnllw8KgRFQG/BQgkj5il9Mw5c+bzW4nDQtATe+yjT5S5epQxwnfjAqG3Y0pt43HkvDagqvIyrlRBIuIqBQs8mS+wUytX7cWK1YsU3lGMTEx8Pb2gbW1FR54cCS6du1WrMeQqTc3NzckJCTA1dW1wveZyl98SjaGfnVB18ft0S6eeK2/Hw91OcrMSEdExGX4+4fA1s6+Qo/t9tPJeGGuJHVr3truauaCKSPYcoaISsekpt5KqnOXLmjWvMVtl0vgQ+YhV5rdLrrZ7LZVTQe8cheb3VZnnes7Y+bDNVRtJRlNWnskCXY24WqUkPlmRFRSZh0oOTk5qxOZL+nhJsvK9ZvdMqel+uvaQBMsSW0lCZZWHEhAPT87PN6NldWJqGSYo0Rm63REOr7ddLPZ7WcPBMHPjUVITUW3hs6YMiJQd/6LtZHYejK5SveJiKofBkpktr5YE4ncGxl6z/byRsc6bHZravo1c8Vzvb3VtmRjykq485GFlxAhItLHQInMkvQH0065Bbhb48nunJIxVRIE92nioraTM3Lxwm83Gx0TERWFgRKZZQL3F2tulgIY29cXdjb8VzBV0shYGhpLAVFxKSYTE/64imyWDSCiYuCnA5mdPw8l4sQ1zfRL40B79G/Bsg6mTgqIfj2qBjycNI2Nd5xNUVOvRERFYaBEZiUjKxcz1t/8gBx/j68acSDTF+Rhq6qtW9941/tteyyW7ouv6t0iIiPHQImqBamLevHiRfz33wFk3OjnFxcbi0OHDiI8PLzYj/P7jjiEx2tqJnWt78QEbjPTtpYj3rrXX3f+g+XhOHAxtUr3iYiMm1nXUaLq4+uvZuDEiROqAXGtWrVx/NgxfPnlNNSuXRs9evYqVmNiqcD94+ab5QDG383CkuZoeHsPnA7PUL39snOAl+eFYcHztRDoztIQRFTGQOnChfOIjCzevL58mPn68oOoPEiLlStXriAzIwO169SBre3NHmTXrl3F9YibicnCPyDAIHA4f/48rKwsUbNmqMHtYqKjkZySfNvl5eXs2TNwdXUr8+tA2sv8++82/Dp7LhwdHdVlq1atwEMjR6Jfv7uL/Tjf/RODpPRctT24jRvq+VdsKw0yXhMH+OFcZAb2nE9FTHIOXpp7BXOfCWUTZCIqW6D096aN2Lx5s+58ZmYGsrM10xjy4Z2ZmanbfmbMcwyUysGZM6cx9fPPYO9gDxtrG8TGxuC1iW+gYcNG6vpNmzZiy+Z/DIKdbt176AKlr76agZMnTyAnOwcdOnTA4088pS7PycnB559PxpNPjUZF+f3339C2TTsMGDioTI8TGXkdzs4uuiBJXL8eAV/f4vdjuxKbifm7YtW2vY0FXujjU6Z9oupNqq9PfygII765iLC4LJXc//bia5j6YBDbnBBR6QMl+VDVfrDKt/w3Xn8VIx54CJ07d4G9vT0S4uOxavVKHD1yGF26dC3JQ1MBZs6Yjg4dOuKJJzUBzuLFi/Dl9Gn4ZtZ3sLLSrOBp1qw5xr487rb7yojRgf378P0PP6tRqdFPP4EHHnwIDg6OWLF8GRo3aYp69eoXeexzc3MRERGB5KQkhNaqpRvROnnihBq9cnd31wVfhw8dQuMmTVRwI7e/ejUM/x04oAK9Ro0aF5qDJLlGiYkJqFEjGM7OmtYy0dFROHbsKHJzc9TjaCUnJ6sRTuQBDRs1VL9TYWasi1LTLNqmt6zATe5O1mol3EPfXkJqZq7qCVc/IAbP9NQUqCQiKlMy97ZtW9TIRe/efVSQJNzc3fHww6Nga2unRkKobJKSknDt2jV06Xoz6OzWrTuioiJx5vTN45uWlqaCicuXLulG+ERKSgpcXFxVYOPk5AR7ewekpqYhLOyKmsp64IGHitwHmdp77dXxeP+9t/Hrrz9j/Lix6nnEd9/Nwonjx3S3lRHFjz/+AAkJ8Thy+DCio6Nx7NgxrF69Epv/+bvA50hISMC777yF9959C3Nm/4pnRj+JNWv+VNddunQJe3bvUr+jPI72lJ6ejl07d6rthITEQn+Hw1fSsOaw5jaeTlZ4gv2+6AaZfpU2JxY3Fj7OXB+Fv48n8fgQUdmTuaOiouDh7pHvdRI4RUVGFjqCQEWTqSZra2s1MoT6DW4c92hdANOwUSMEBQbh4oULWPy/ReoySysrvPzyODRo0BABgYFIS0vFyhXLkZ6RDkdHBzX688XUzzB69BgV2Jw/dw41gm+O4NxKRq9q1a6NMWOeUyNYMsITfWMfCnNP/wHYs3d3sabe/pg3V40ofTPrexXUHTlyGB99+D6aNG6KNm3awtLCErNmfY133n1fd58nn3gUIx9+BC1atCz0seVx9evlSDsLZ3vNSByR6NXYBS/d6YMZ66PU+YkLr2HeszVRnzlsRFSWQKlu3XqY+9tsNG/RAvVvfIjLh9K2rVvUku1Rjz7OA1xGEpj07XcXfvnlZzU6JEHTypUrVDCRnaOZR+rVu486aae+Zs/+BVM/n6Km22xsbPD2O++pQEkCqDffeherV61UfzsnZye89OJz8PcPUPk+73/wEYKDQwyeX/Khzp07iwmvvqab5vP29lGn0pCRIwnMtGSKzs7ODgcPHsTDj4zSTenJVGJISE0cPnwIITVroiz+OZGMfRc0y79DvW0xrH3+wT2Zt6d7eOF0RIYaeZRpuBd/k5VwofBw4sJgInNX6neB7t17qGmXN9+YCC8vL7i4uKi8JZkSeezxJ1CjRo3y3VMz9fjjT6pEbcn7srCwUCNB77zzJnx8bg9WJJjp338g/vpztcopCgwMVPd98aWX1fUyjbd58z+YPOVzzP1tDu6+uz/uHz4CS5csxpq//sLoZ8YYPJ4EZ0Km7YrnRofZAkjekkyVaclIlQRKMuqln6gtHJ0ckZpatvo20qJi+tqbo0mv3OWjkniJbiX/Wx8NDcCl6Ewcv5auErxf+eMqfnwihK8ZIjNX6kDJ0tISzz3/IgYMGITjx48hMTERXt7eaN2qNTw8Pct3L82YFFfs0+dOddKuPJTApXHjJrpgRj+Qkek34ebmZvA4Mtr37ayv8eRTT6vgJCUlGXXq1lXXeXt74+Kli7c9t4w2yW2PHDmCjh3v0F0ugY0kT8t0XWLSzfygK5evGNzf2tpGJYJrSeK4/vSZlgTVp06dVNNsmsdPw6WLF3HXXfegLKTq8vkozUrMVjUd0LuxpjEqUX4cbC0x85EaGPHNBVUyYO/5VExefR3v6BWoJCLzU+ZxZZkaKev0CBVsw/p1KnlbAiOpTC0jMi+8OFYFMOKD999F23bt1FTV9YgIrFixDPfc0/+2UaC//lqNkJAQNG3aTJ1v0rQZ/ly9CvZ2dli9ehX63Nn3tueWqbuHRj6iAqyI8HD4+Ppi3949aN26Dbp2667yg5YvWwoHewdkZWdjzV+rDe4fFBSEPXt2q1wpCaoKylkbOmy4ypuytbFVZQ3WrV+rRsxktV9ppWTk4OuNmpwTMeFuXy77piIFuNtg5sM18NiPl5GVk4cFu+LQwN8OwztwypbIXFnkyVBDGciSbqnTI6UBvH280bx5S7i6mk+TURlJk9Ebyb+piN9b/jwSLB08dBDubu7o3acP6tTRjARpR5TWrvkL58+fg7OLixrRa9+ho0FQIKNSskJNpu0cHBzUZTLS89dff6rpU0kKlyk7GSXMj+Sc7dq5AympqWjXrj26du2mLpcVdjLNJ6NBMip1Z99++PWXn9VUnySNy6q9ZcuW4GpYmDr/7HMvFPh7nj51Sk0LSnkAKUEgI5Xa1ZSSJyXB3Nixr+huP336Fxg8+D5V2DQ/32yMwqxNmqTzvk1dMH0kp4KrUmZGOiIiLsPfPwS2dsZf6HPZvni8vUTTGkemaxc+H4oGAca/30RkZIGSFBSURGE7O3u4urogNjZWTbc8/fRoVTrAHFR0oEQlF5WYhbunnkNaVp5qgLrylTqo6X2zmjlVvuoWKIlPVkZg3s44tV3Pz04ld9vbsD0mkbkp9dTb/v371EiGVIlu27adukyKGq5d+xe+/fYbtXKJuUpUFb7eGK2CJDGigweDJCoV6QW490KqWg135nqGWhjwxkDmKxGZm1J/PTp86CD6DxioC5K0OS0DB96rclGOnzgOYyeB3dYtm/HLzz9h3u9zVf0eqt7OXs9QSdzC2c4Sz/ZmlWUqHTsbS3w2IhC21ppp7N93xGHbqWQeTiIzU+pASSbsrCzzL9wn+TF5eTdXOxWHlBWQlVj6J23vuIog+TNjX3oB/x38T/Wks7C0wGdTPlXVp6n6mrY2Erk3JpOf6uHFOjhU5srdshBA663F1xCTfLP6PRGZvlJPvTVr3hxffzVTrZ5q0qSJruChNGmVcgFSOqA4pKK0rOSSwOVWw4c/gBEPPIiKIMUNP/r4E1UDSqtmSE1MmzYV9947GJ6eNy8vjgvnz6taUlpOzs7w8/NTwV7YFcNl86J2nTrqp/RCy0jPMLhOVpfJY0nek6rKrcfewUHVR5JjLUvobyUrEKUwZUREOFJTDOsQeXp5qaRq6ZMWef264fGws1U91sT589JDzTB1LahGDbXSLjIyUvVw0yeta+Q4StmA8GuaBFgtK2srXcPeS5cuqua8+gICA1SpAanBJQsC9ElyugSxkowuCeEGLCxQu7YmkVtasmRmZOLo1TTsOhgBWQ/o6eWDRzp7Ij4+HrExMbfVaJLSB5KMrm3Hoq9maKiqSSV1p9LT0gyukxIYkpMmr1epPq/Pzt4OQUGapHH9wppaUgFdXnfXr19HSrLhyISHh4eaqpbaUbLCUJ+1jY1asShk5WNuzq3HMFAl6ctrRV4z+lxcXdUKwvyOoXw50CbDX7lyGVmZWQbX+/n7wcnJGXFxcYiLjb3lGDrB399fjcpeuXz5tt9VEvJlcYCUq0hPS1c1tqyt5W8VBm8fX5XPJ/l90VE3VyYK6QkYGBikFhtIxflbBYeEqJFrqROWeqPOl+4Yenqq4yhfsq5HGL6+bWxtdAVVpUdgnjaavuX1LR0HkhINW+L0q++KrfWdsP1kLJKjo/HGL9GY2N8XFrBQhVxDQzWv78uXLyM7y/AYSi9EqREmx0+Ooz6+R1T+e4Q+Xz8/tRqX7xHG8h5xk7ePT6W9R2g/iyskUJLVT1J0UnqAubq6qV8qJiZafYA/9fQzBgFIYXlOK1cuV1W8Q4JDYGllOMDlWYH1mORNUbvEXv9NTSQmJpU4UJIikPLH0e/JJo1q5Z/71Vdvb1i7ZOkK9VOCzdOnTxlc99LYV9Sx3bHjX/z04w8G17Vo2RLvvvuBelHn97i//Pqb+iCX1Wf79u01uO7Rx57AoEH3qorXshxfn/wzTP1iutqWZsf6PePE9C+/Uh/Wi/+3UAXD+u4bMlT1+Dt37hzee/dtg+vkOP740y9q++OPPlTVvvV98OHHqmSB9HZbtnSJwXXSR1ACbgksbv1dJRhcuGiJrs2KapArv8eN67u1eR72Ni2wfusWzJmteX4tmS5+48231YrB/I7h3N/nqw+3n376HocOHjS47qmnR6tCnQcO7FcNi/VJhfpPJ2uOa36P+/U336nyBwvmz8PWrVvy/VIgKwg//ugDg+vkzUbau6jj9f476s1D3yefTlEta1auWqEqr+u766678fToMeoN8NZ9kuDq93kL1PbUqZ/dFtC//vqbaNe+A/75exPmzZtrcN0dd3TChFcnIjEhId/fdcHCxepN8LtvZ6k+hDVqBGH8+JcxY8Y0DBw4WJWjkB5+ks+or0mTpvjwo0nqfSS/x/3hh59VsPr73DnYuXOHwXUjRz6CIUOH4fixY5g8+ZPbgtQZM75W2++8/aaq1aXv88+nqTfM5cuWYO3aNQbXSQuej4c9iuGTTsLz0jeIuQS8tllznbzv/Tpbc2ymTJ6k3pz1SWX8Vq1aY/36dVi0SHOstfgeUTXvEVrjJ7yGTp06q/9FvkdU/XuEvmeffb7S3iO0n8UVWh5AvkWdOHFMfZP19vJWbwrFTeKWF318fBwef+IpGAP5g+zftxffff+T+ifLj0THctKSDy358Dv4338cUarCEaVNR2J1vbpqetlg9kut4ebqwm+LRjWilIvsbMtqN6IkXzzkjXf9wSi8PueQuszO2gKTRwSiprcDR5Sq0aizPo4oGdt7hPGOKJU5UCqL7dv/xe7duzBu3IRyeTwpVXDp8u1TKfqkuav+yI+WjJJ89+03aqRBCioWZOGC+QbfDOVFsGr1GpYHqEKZ2bkYMO08rsZp/ol/eDwYnevn3+SXqkZ1LA+Qn49WRKgilKJhgB3mPxcKW6lBQUQmq8yVuctScLJ9+w74c/VKLF++FG3btofDjQKDWg6Ojrf1ACuMfCuRbxiFya+o4r//bsUP33+L5194qdAgSciw3cBB9942okRV54+dcbogqVM9JwZJVGEksXvP+RScj8zEyfAMzNwQhQl3+/GIE5mwMgVKZS04KavjXFxcVYNWOZU1mVvyd0pq+7/b8NXMGRjz7PPo0aNnkbeX0Sj9Ealbc3mociWn5+D7fzQJ71KMXGrfEFVkP7jPRgThgVkXIDNEv26NRZd6zuhYt7iNo4mouqnSgpM7d2xXyasy9Sb5A7KCRF9FV7resWM7Zs78UgVJPXv2qtDnooqx4kACEtM0pSgGtnRDQ7aZoArWKNAer/Tzxed/aVY9vvm/a1j6Ui24O5V5gJ6IjJB1RRSc/O/AAVVwsnPnLoU+hiTnde/RE527dEVlk4TDL6d/oVZcHD16RJ20Bg4YpBLNyLhJet2C3TeXXD/RreJWSRLpG9XZUxWf3HUuFdcTs/H+8ghMfyiIjZeJTJB1VRaclPovZ8+eQVVwcnJUI0n50a+HRMZL2ktIrohoW8tRFQckqgyWlhb45P5A3DfjAhLScrDhaBKW7U/AkLbu/AMQmZgqLTgpdWd+m/MrFi1cgHbt28PeXtPZXkuKgVVU0CJLGnv16l0hj02VY/6NhqXigY4ePOxUqfzcbPDhUH+M/f2qOv/Jqgi0CXVkb0EiE1OlBSfXr1urah0sXDhfnSqzMjdVb5GJWdh0XFOrxcvZCr0bcxSQKl+fJq4Y1i4Fi/fGIy0zD68tvIrfx4TCxkrTH46Iqr8qLTip6elmWBBKn5OTkxr5MWZSHkAK0snvX9HJ53TTrI1R+GaTZrXbMz298FJfrnYzZqZSRyk/qZm5uP+rC7gYrZkGHt3TC2P5eiQyGWVepiFtLbR9qEpKgiBjD4TI+GTl5OF/ezUVeqXrzfAOnHajquNoa4kpIwIx8tuLyM4Fftwcg871nFXeHBGZeaAk5QCkZ4pMuUk5cX2NGjfWNQgtiDTxlJVvBZHpO2kdQKTvnxNJiEzU1K/q2cgF/m63V1onqkxNazjghTt98OW6KLXQ5fVFV7H0pdpwdch/wQsRmUGgJFNNb77xGqKjo9XUk+UtK+AeeWRUkYHSxo0bbmsUqY85SpQfJnGTMXqimxe2n05RqzHD47Px4fIIfP5AIEsGEJlroLRxw3p4eXnjs8+nqVyi0pB2IIPuHWxwWXpaGnbu2olNGzfg3sH3lXb3yESdi8zAnvOpajvU2xYdanN6g4yDlaUFPh0eiCEzziMxPRdrDieiWwNnDGrtVtW7RkRlUOpujolJiejQsWOpgyRtgUoHBweDkySC33NPfwQFBeHQoYOlfmwyTQtvNCQVIzp4qHo2RMYiwN0G7913s/ejlAyISjTsuE5EZhIohYbWwuVLl1BRpGJ2RHh4hT0+VT8pGbmqZYmwt7HA4Db8pk7G567mrhjYUrMCNik9F5+uvl7Vu0RElTX1Fh4ejvh4zTd6Hx8fLJh/EHPnzkHzZi1ga2drcNuAgEC4uxdepTYzMxNZWZoltVo5Obm4cOE8tm7djGcLqJxN5unPQwlIztAsGujf0o2JsmS0Jg7ww79nUhCXkoN1R5LUAgRZeEBEJh4oLVu6WFXe1rd82VJ1upUEOX3u7FvE4y3JN5lbWqBI1ey27dqXZPfIhEm5L4Mk7g5sFUHGy8PJGq/398PERdfU+Y9XRKB9bUc42XEVHJFJF5yUqttyKg4rKyt1KkxsbAxiY2MNLrO2soa3j49qX1IdsOBk5ThwMRWPfK+Z6m0R7IA/ngutpGem8mDKBScLIm+to3+9gh1nNEV1H+7kgTcG+lf1bhFRRY4oFSf4KWkekpyIirJAL4mbfd2oOpCR8XcH+2Pwl+eRnpWHeTvj0L+FG5qHGPa0JCITCpSkXYkUiZRK3CmpqWq7IHKbkhSLTE9PR06Opoiglo2NLWxtDXOfyPzEJGdj3dFEte3hZIV+zZjrQdVDsKctXujjg6lrIlUhyveWhWPRC7XYC47IVAOlv/5cjU2bNmDMs8/jzOnTarsgcpvevfsU+ZiL/7cIq1evRFKSpsGpPhacJLFkbzyyb8z4DmnrDjubUi/WJKp0j3T2xOqDCTgZnoHTERmYvS0GT/dgxwEis2mKWxb79+/DjC+nYdSjjyMkOASW0rhLj6enp9FPzTFHqWLl5Oah72dnEZGQDQsLYO2EOqjhyVHG6sYcc5T0Hbuahge+uYjcPMDO2gLLxtZGTW++jonMoiluWUgdpp49e6FPnzurcjfIiG05mayCJNGtvjODJKqWmgQ5qJGlOf/GIiM7Dx8sD8fPT4awvQmRqeYoFUdxcpR8/fxw4eKFkuwCmRkmcZOpeL6PDzYcTcK1+CzsPpeqiqcObsMyF0QmmaNUHMXJUWrfvgP+XL0Sy5cvRdu27eFgbzgk7+DoCEdH9vIyV5eiM7H9xtLqGh426Fy/9O1yiKqak52lWgU3ZvYVdf6zPyPRtYEzvJyrdGCfiIw5Ryk7OxuffzYZ+/btzff66pDMzRylivPZn9fVVIUYf7ev6s5O1ZO55yjpe3XBVfx1SLOKc0BLV0wZEVTVu0REhajSrzI7d2zHqVMnMW7cBASrZG7DGk2urpp+SWR+0jJzsWx/vNq2tbbAfezrRiZiYn8//Hs6GYlpuVh9MBEDW7mhS/3qUWCXyByVKVDKzc1VwY7kGbVr1x4NGjTEtWtXYWVlDT8/vyLvHxMTg+49eqJzl66oyhEhCdZkdEsa/QYE3Oz8TVVn7eFE9UEi7mrmqlpCEJkCbxdrvHqPH95Zomn6/eHyCCx/uTYcbVn2gsgYlfrTR1qZfPDBuwi/Jr2MLODv568CpYyMDHz37XRMnvJ5kSs6AgIDcfbsGVSVRQsX4J/Nf6NWaC31+xw5clgFbqNHj+FqlCo2X68S94MdPap0X4jKm4yQrvwvAXvPp+JqXBZmbYzChHuK/nJJRNUoUNq2bSsyMzLw9Tff4Zeff9RdXqtWbVXvRkZpGjZsVOhj1K/fAL/N+VUFLO3at4e9vWFpf+n35uJScVWYa9epg6HD7te1ZTl27BjefedN9OjRUwV9VDWOXEnDsavpartxoD2aBZt3TguZHvkS+f5gf9w38wIys/NULt49LVzROIjtTYhMJlC6cP4cunTpBjs7O/mvN7jO188fEeHhRQZK69etRUREBBYunK9OlZ3M3bZtO4PzXl6aZGGZhiMjGU26w4Oje2SSQn3sMKaXN2auj1KFKN9bGoH5z4XC2qrwkXgiqiaBkrWNDZJTkm+7XBbRXblyWY3KFGXAwIHo2atXgdc7OVX8cvBr167h+LGjSE5OxtatW3DX3fegSZOmBd4+KytLnbRSU1MrfB/NSXxKNtYc1qwIcrW3xN3NmdBPpuvxrl5YcygRZ65n4Pi1dPy+IxaPdeXqTiKTCJQkeXvq51PQrVt3WEDzDSgzMxPz/5iH2JhYNG7cpMjHcHJyVqfy8t9/B4osiNmzV2/dVJtITEzA6dOnEB8fr7a9imiZsnTJYixatEB3Xj9oorJbtj9BTUUIKcbnwARXMmGyovP9+/zx8PeXVNPcrzdEoU8TF1agJzKVOkqLFy9S+UU2NjZwcHBASkoKLC0tMWHCRLRq3Trf+8TFxSE+Pg7u7poEXdkuiNzGw6P4ibyrV63E5cuXCr3N06PHqP3Nz/nz5/H6xAkYN/5VdOx4R7FGlGTVnKyUS0hIYDmDMsrNzcM9X5zDlVjN8f1zXG01PUHVH+soFe7jFRG6Kecu9Z3w3WPBnHImqu4jShJfDRs2HJ06dcaBAweQlJgITy9PdOhwB9zd3dXqN5W/lE9ekozISP6R0B+dKWuO0oCBg1AWtWvXhp+/v0pELyhQkiBLP9BiPlP5kSrc2iCpU10nBklkNl7u54O/jyfhemI2/j2dogpS9m/pVtW7RURlCZRWrFiGnj17IzAwSJ30bdq0Uf3Mr4XJ4PuG4O57+uuCKNkuSH6BVnlJS0tDdnYWXFxu5sDExsYgOioKfn7+Ffa8VPwkbiJz4Wxvhbfv9ceLc8PU+cmrr6NzPSe4s34YUZUrdYWzjPQMfPzRB7clM2/cuAE//fg9atSoUWDwIxW35aeMxsh2ficpCyB5QxVFRrzefusN/PTjD6qHneRWvfH6a6hduw569iw4wZwqxtW4TGw9pVkc4O9mjW4NWKmYzEuvxi4qP0nEpuRg6prIqt4lIipLoDTs/uHw8vbG5E8nqSRusXHDevz80w+Y8OrEYtUh+nP1Kvzw/XcqaNEXFxuLjz56H5v/+bvC/kgyPTh5ylTUDA1FeEQ4LCwtMHr0s/h40qcVOpJF+Vu0O14ls4rhHTy4RJrM0luD/OBsZ6lb2HD4clpV7xKR2StTMrcESDKq5OjkhFatWmH2r7/g1ddeR+vWbYp1/4sXLuDzqVNgbWWFl18Zr4pV7t61E99++w1CQmripbEvw9vbx6j/SGyKW3aZ2bnoNfks4lJyYG0FbJxYDz4ubFliSpjMXXzzdsTik1XX1XbzYHvMGxMKS0vWViKqloGSkKm3d995C2FhVzBx4psFrnYrLFdIKntLDaNmzZqrNiIjRjyocplkBZ2xY6BUdpLEqs3NkLpJUx9kN3VTw0Cp+LJy8jB05nmci9SM1H96fyAGtWZiN1FVKdHX9pUrV+DAgX35Xic5RStWLlMnMWjQ4GKNLElZASnyuG/fPlUHqUWLlirBuzoESVQ+NhxN0m0PbMUCk2TebKws8PoAPzz9yxV1fvq6SPRu4gKnG1NyRGTEgZKvry/q1KlbrNu66q0mK0hubi5WrliO+fPnqRV0vXr1VtNuE8a/jLFjx6F+gwYl2T2qhqS45D8nNIGS5GbcUbfiq7ETGbtO9ZzRs5Ez/jmRjMjEbPy0JRpj+/pW9W4RmaUyT72VxYrly7B8+VI8+9wLaN++gy7vae7cOare0ktjX0Hnzl1gzDj1VjbbTiVjzGzNN+cBLV0xZQSn3UwRp95K7lJ0JgZ9eQ7ZOZoK3qteqc2K3URVoErHckNq1sS0aTN0QZKwtbXFk08+jddffwspybf3kiPTsuHYzWm3O5tqlkYTEVDT2xaPdPbUjbyyXABRNZh6kz5ne/bsxtChw3DlyhW1XRC5TTu9ACg/rVoVnPhd0qRwqn6yc/JUIrdwsLFA53qsnUSkb0xPb6w8kICY5ByVy7fnfAra1+b0NJHRBko1a4aqvCJvHx9YWlqp7YLIbYgKs+9iqioJILo2cGYDXKJ8Kna/3M8X7ywJV+c/XXUdi1+sBSuWCyAyzkCpTdu26iSk5pF2uywuXryo+r1dunTxtqm2wYOHqDIBZPqr3TjtRpS/wa3dsGBXHI5dTcfpiAws2RuvirISkRnkKCUkJOD9996Gh7s7aobURJMmTTFk6P3w8PSEk5NTsQtXUvWTm5uHTTfykyRRtXtDTrsR5UeKTU4c4Kc7P2N9FBLSNCOxRGSkOUrFUZwcpQMH9qtWJ0+PHoOFC+YjD3kYNOhe3HXX3Xj11XFISr454kCm5eDlNEQlZavtTvWc4GRnVdW7RGS02oQ6qmKsaw4nIj41B99uila1lojISHOUiqM4OUrRUVGq15qws7dTPd60K99atWyNM6dPq1EmMvFptxuNQImoYOPu9lU1x9Kz8jB/ZyyGt3dHbV/2pSQy2hyl8iBBl5WlZiTBz88fO7ZvV5dJVW5pieLnz29MpkhKd228Me1mbQn0bMxAiagoge42eLKbF77ZFI3sXOCzP6/ju8dDeOCIqkuO0ratW3Dq1MkS3cfC0kKdRJs2bZGUlIhXXn4R48e9rB6rffuO5bV7ZEQkKfVafJba7lDHCW4OnHYjKo7Hu3nB303z/Xbb6RRsPclac0QVrdxatEsz2/oZDVTOUXENH/6Ablum2yZPmaoCLqnO3alzF3h5eZXX7pER4Wo3otJxsLXEhHv8MGH+VXV+yp/X0bGuk1oQQUQVw6i6LLq5uWHAwEEYMnQY/P39q3p3qIKm3bSBkgwm9ua0G1GJ3NXMBa1DHdT2xehM/LFTk9tJRCYeKKWnpyMlJdngJCNLZFrOXM/ApRjN37VNLUd4OpfboCaRWbCwsFAr3ixuDCLJCriYZM0KUiIqf+X2KTXq0cdhY2NT4vst/t8irF69EklJSflOzY144MFy2kMyBlztRlR2TYIcMKSNG5bsS0ByRi6+2hCF9+8L4KElMuZAydm55AUD9+/fh5Url6sgKyQ4BJZWhgNcnp6ahpBkmoFSH5YFICq1l/r6Yu2RJKRk5GLx3niM6OCBRoH2PKJExhIoSduR3bt25XudLO93cXFBs2bNcU//AbCzy7/Wx+VLl9CzZy/06XNnaXeDqpELURlq6k20DHGAn1vJRyCJSMPbxRrP9vLG1DWRyMsDJq++jtlPh6ipOSIyghwlWd2WmJiAzMwMtGjZEl26dEXDhg0RGxsDa2trVZzyr79WY+rnUwp8DF8/P8TFx5d2F6ia4Wo3ovI1spMnQrw0Xzj2XUjFer0RWyKq4kBJkq/9/QPw5YyvMWrUY7hvyFDVimTa9JmIjY1VI0lfTJuBM2dO4/z58/k+Rvv2HRAdFYnly5ciLCwMMdHRBqfU1FRUlrS0NLUP+eVKUflgoERUvqQswGv9bxbm/fyv60jPKl73BCKq4Kk3KQgpvdysrAyLBXp4eKBe/fo4d+4sOna8A42bNEX4tauoXbv2bY8hQ8QuLq6Y+9scdaqqZG6pBj5l8ieqFtTIhx/BkCHDKvw5zU1YbCaOX0tX202C7BHkYVvVu0RkEno0dFb9EnecSUF4fDZ+3RqDZ3sX3UKKiCo4UJI8pJMnT6gmtvoyMjJw4fw5dO/WXZ1PTk6Cl5d3vo+xc8d2FXCNGzcBwSqZ2zDocnV1RWVYunQxHB0d4e7uXinPZ460LUvEnU3ZsoSovMgXzon9/TBk5nnk5AI/b4nBfW3d4c8cQKKqnXqTJOxDB//DpI8/xN+bNmLPnt34668/8cbrr6kgqnmLlrh48SLS09LVCFN+YmJi0L1HT3Tu0hUhNWuiRo0aBqfKCJQk2Nuwfh2efe75Cn8uc8ayAEQVp66fHR7s6KG207LyMG1NJA83UVUHSkFBNfDp5M9V65G5c3/D559NxvJlS9C4cRN8/PGnaqWbn58vPvr4k9um57QCAgNVLlJVSU5Oxowvp+HZ515QU4BUMa4nZOHg5TS1Xc/PDqE+7HhOVN6e6+MDd0fNe+2fhxJx/Krmf46IqrCOUkhICF597XVda4pbl6U6ODgWev/69Rvgtzm/YtHCBWjXvj3s7TVl+fVrM0mZgeKSoCstXZMHU5CgoCDdfs6a9bVKKG/ZslWxnyMrK0udtCoz4by64rQbUcWT5tJSLuDT1dfV+enrovDjEyE89ETGUnCyNLU71q9bi4iICCxcOF+dyprMLbWdjp84Xuhtpk6drka7Nm/+B2dOn8bw+0eo1W7apO7EhESEh4cjICD/KrdLlyxWz6OlHzRR/jYwP4moUgzv4I7ftsfialyWSu7edS4FHes48egTlYFFngwFlVJUVBQWLPgDx48dRXx8PLx9fNRKtyFDhhY5miQ0Pd1SCrzeyckJTk4lr/hdHKtXr8K6dWsMLosID4ezs4tqyPvp5M+KNaKUmJiogqqEhIRKSz6vTmKTs9H9kzPIzQNqetniz/G1WRDPDGVmpCMi4jL8/UNga8fq0RVp1X8JeH3RNbXdtIY9FjwXyv85oqoYUUqIj8fE18bD19cX/fsPhIurK6IiI7Fp0wYcO3oEH0+arJK6CyNBUEUFQkUZMGCgOul78olH0X/AwELLA0g/O/2edtnZbEZZmE3Hk1SQpF3txqrBRBXrnhau+GVrDE5HZOBoWLoa0e3blF/iiCo9UPr7702q+vY7775vEBD1u+tujB83FidOHEeTJk1LvWNkeqvd+jZjWQCiimZlaYGX+/nguTmalIKZ66PQq5ELrK3Y2oSoUle9ybRbq9Ztbhs1kuTr+g0aqtGl6iYgIJCr38pRQloOdp/TTK0GedigMRt2ElWKbg2c0SZUszjmQlQmlu1nqyiiSg+UvL29cfTI4dsul1VgkiTt5Z1/kUlj9vGkT3HnnX2rejdMxj/Hk5B9o5tCnyacdiOqLDLF/cpdvrrzszZFs7UJUWVPvUnByRUrluP9995Bp06dNTlKUVHYsGEd3NzcVD2l0tY2kma7MrrDfJbqjavdiKpOq5qO6NnIGf+cSEZkYjbm7YjDk929+CchqqwRJQ9PT0z57HO10kuW9k/74nP8uXoV2rRui/fe/6DAIpO3+ubrmbh44YLalga6z4x+Ei++8Bw+/XRSaXeNjEBKRo5anix8Xa3RItiwRhYRVbyX+/lCW7nlp83RajqciCqxjpK/fwDGjX+1wIKTRZEmtNLGJLRWLV2Nol69++Cuu+7GJ5M+wuHDh9C8eYuy7CJVkS0nk5GZrVnu1ruxCywtmUhKVBWtTQa1csOKAwlITM9VfeDG6U3JEVEFjijdqjTTZOfPnUOtWrXVdk5ODo4ePYJBgwar9ijtO3RUveLIBHq7sQkuUZV5oY8PbG6sePt9e6xqKUREFTSiJCM+0vy2OIYOHYZ27TsUehsHR0ecPn1KbR8/fgyubm7w8fFR51OSk1XARNVPWmYutp1KVtseTlZoE1p08VEiqhiBHjaqYa5U7M7IzsO3f0fj/fvy7zxARGUMlKRukrT5KA6p0l2Udm3bYc7sX/DRh+/j8uVL6Nfvbl0RR6nDNOz+ESXZPTIS288kqw7moldj1m8hqmqje3phyb54pGTkYum+eDzaxRO12JyaqPwDpTZt26pTeZGE8E8+nYLt2/9F+w4d0KePZmn+xYsXVLVvqfpN1bzIJKfdiKqch5M1Hu/qia83RiMnF/hqQxSmPcQRe6JKbYpbWjJKJSd9devWUyeqfjKzc7H5hGbazdXeEu1rsyEnkTEY1cUL83fFISY5B+uOJOFoWBqa1uBqVCKjDpRioqPVqreCeHl5VcvCleZs19lUJGdopmd7NHKBrTVXuxEZAyc7S4zp5Y1JK6+r89PXRuLnp2pW9W4RGb0qDZQ2btyARYsWFHj98OEPYMQDD1bqPlHZrD+aqNvmajci4zKsnQd++zcWV2KzsOtcKnacSUanelXTmJyouqjSQGnI0GEYdO9gg8vS09Kwc9dObNq4AfcOvq/K9o1KLisnD38f10y7OdpaonM9TrsRGRMZ4X3hTh9MXHhNnf9yXRQ61nFinTOiiq6jlJGRgaysktfmsLGxgYODg8FJErzvuac/goKCcOjQwfLYPaok+y6k6ir/dm/oDDubcivTRUTl5J7mrmgQYKe2j11Nx3q9xRdEdLsyfZIdOLAfL734PEY+NAJbNv+jLtu7ZzeWLV2CsvL09EJEeHiZH4cqzwZOuxEZPamS/0q/myuKZ6yPVKPBRFTOgVJYWBimT/sCffvdhTZtbpYMaNmqNdZvWIeEhIQiHyMzMxMpKckGp8TERDWStHXrZgQGBpZ296iS5eTmYeMxzTdTO2sLdKnPvAciY9WlvhPa1dIUgr0ck4Vl++KrepeITC9HadfOHejVuzcGDBiIK5cvGUynhQSH4Pixo7ijU+dCH0NGnvJL5pZ2KL169Ubbdu1Lu3tUyY5cSVPLjrVvwrLChoiMk7zHvnKXLx76VtMmatamaAxs5QYHW/7fEpVboBQXF4vAoCDNmVv6vOXdqK5dlDv79r2tgKW1lbWq6u3szBGJ6mTn2RTdtpQFICLj1iLEAX2auKiR4KikbPy+IxZP92A5FqJblfrrQ1CNYBw7dkxtW+BmoHT9+nUcO3oEobVqFSsPSVtcUnuS+zFIqn52n0/VbcsqGiIyfi/19YHljbfvn7fEID5VMypMROUwotSjR0+sWrkcM76chpjYGDWoJMUj16z5Ey1atkJwcEiJ8p2ioiKRk2P4Tyo5SoGBN0atyKib4B68lKa2gz1tVBNOIjJ+dXztMLiNG5buS0BSei5+2hKNCXf7VfVuEZlGoOTo6IiPJ03G73Pn4OKFCzh29Cjc3d3Rs2dvPPjQyGI9hkzPTf50Ev777wCsrKxgaWk4wDVs2HAMu394aXeRKsl/l1J1q2Y6cDSJqFp5rrcPVh9MRGZ2Hv7YEYeHO3nC341fdojKpeCktBgZ+/I4tS2jQRLslMSuXTtw9WoYpn/5FYKDg1WCIVU/u8/pTbvV5bQbUXUS4G6Dh+7wwOxtscjIzsOsjVH4cChXHBNpldsSh5IGSSImOgYd7+iEkJAQBknV2K5zNxO529fWLDkmourj6e5ecL6xUnX5gQRcis6s6l0iqp4jSkuXLMaePbuLdduhQ4ehXfsOhd6mRnAwTp06iaoyc8Z0VTRTX+MmTfHaa69X2T5VN4lpOTh+NV1t1/e3g5dzlXbFIaJScHeyxmNdPfH1xmjk5AI/bI7GpGEcVSISJfpUq1kzFLm5ms7wRZEl/kWpV68+5v0+F4sWLkC79u1hb+9gcL2sfnNxqbil5qlpaWjfvgMefuRR3WXW1vygL2nbktwbRX2Zn0RUfT3c2RO/bY9FYlouVv2XgNE9vFHT27aqd4uoypUoKpCaR7fWPSqLNX/9iUuXLqrTwoXzb7t++PAHMOKBB1GRbGxt4erqWqHPYcp26dVP6lCH025E1ZWLvRUe7eKFrzZEcVSJSE+VDp8MGDgQPXv1KvB6J6eKTwzeuWM79uzeDTc3VzRr3gL33z9CreijktVPsrKEriUCEVVPIzt5YM6/MRxVIiptoPTXn6tx8OB/uKf/AFy7elVtF0Ru07Jlq0Ifz8nJWZ2qir+fP7p06YL69Rsg/Fo4Zs/+BSeOH8OkT6YUmJyelZWlTlqpqTdXfJkbqeZ79nqG2m4SZA9n+5In9BOR8eCoElEZAyU3Nzf4+wfA0cERFy9eUAFDQYUl5Tb5iYuLQ3x8HNzdPdR52S6I3MbDQ3O74vhi6mc4cuRwobf57vufYG9vr7Yfe/wJ3eW+vn4YP+E1jH3peRUA6jf6vTWhXb8/nX7QZG726K12YzVuItPAUSWiMgRK7h4e6NvvLtSoUQMbN65H585d0OfOviV5CKxft1YFGpJ/JPJrilvaHKVnn3uhyB5z2iApP/J7OTg44OrVqwUGSkOGDsPAQffqzicmJiIgIADmaLdeoMREbiLTGVUa1Zkr4IhKFSht3bJZrVSTgKK0Bt83BHff0x92dnbqvGwXRHub4iprbpGMdqWlpRWa3G1jY6NOWsVp/muqdt0oNGlrbYGWNQ1XLBJRNV8B928sEtO5Ao6oRAUn3dzdcebM6TJNN0nwI4GI/NRuF3QqaaBUEjHR0Zj96y/qp/b8VzO/VNN97dq1r7DnNRVhsZm4Gqd5HbSq6QB7m3KrXUpExjCq1MVTbWvrKhGZK4u8vLwbVXCKdu3aNbz91utISEgo8rbPPvt8iabl0tPTkZNjODpjY2MLW9uKqeMhLVfWrVuL1atWqjwp6TPXvHkLjHr0MZWHVVwy9Sa5W3JMzKnMwOK9cXhvaYSuA/kzPb2repfIiGVmpCMi4jL8/UNga1fw9DcZj6T0HPSdclaNKsmq1lWv1GFdJTJLJZp6CwwMxKxvf8Dly5ewaNFC1AiqgVatW+d7W+ndVhyL/7cIq1evRFJSUqXWUZJVbffc01+dMjIyKnT0yuT7u7ERLpHJjiqxWjeZuxLXUZJkaFlOLwGGp6cXQkNDS/3k+/fvw8qVyzHq0ccREhwCS/naosfTUzP0W9EYJJWMDEJqE7md7CxVaQAiMv1cJRk5DvFitW4yL6UuONm6dZsyP/nlS5fQs2cv9OlzZ5kfiyqP1E6KSc5R221rOcLayoKHn8gcRpX+icbH7AFHZqZMlbmjoqKwYMEfOH7sKOLj41V/t44d78CQIUPhUEAdJX2+fn64cPFCWXaBqgCn3YjMc1RppfSA46gSmZlSL1VKiI/HxNfG42rYFfTvPxBjnn0e3bv1wL/btuKjD98vVvNcaUgbHRWJ5cuXIiwsTK080z+Zc9VrY7bLoH4S25YQmdUKuH+4Ao7MS6lHlP7+exNq1gzFO+++r1aMafW7626MHzcWJ04cR5MmTQt9DAsLC7i4uGLub3PUqSqa4lLJZOfkYe+N/m6eTlao58ckeCJTx1ElMmfWZZl2a9W6jUGQJFxcXFC/QUNERUYCTYpuSHvq1EmMGzdBtUKxvKW/mjktt68ujl9LR3KGZrSwfW1HWFoyP4nI1DFXicxZqafevL29cTSfvmoyXXbm9Gl4eRddVycmJgbde/RE5y5dEVKzpqr4rX9ioGTkbUvqOlXpvhBR5Y4qudprPjIkV+lyTCYPP5mFUgdKslrt1KlTeP+9d1T/tp07d2DlyhWYOHGCKsDYuHERw0kAAgIDdZWxqXpgIjeR+Y4qPcJcJTJDpQ6UPDw9MeWzz9Woz8KF8zHti8/x5+pVaNO6Ld57/wNV0LEoUo/pwoXzWLRwgfoZHh5ucMqvCCVVnYysXBy4qMlPCnC3RrDnzZ53RGT6Hu7EUSUyP6XOUQoLu6Ka0I4b/6quCKEkZ5eEjERFRESoQEtOt2Iyt3E5dCUNGdl5umrcJf17E1H15uqgGVX6hnWVyIyUOlCS0aPaderizhv93ErzoTlg4ED07NWrwOudnJgDY0x2ndUvC8C/DZG5jirNZV0lMiOlDpR8ff1w7erVMj25k5OzOlH1y09i/SQi88RRJTI3pU/m7tUbe/bswsYN61VCdlpamsEpJ0fT4oJMQ0pGDo6Epant2j628HVlfhKRuWKuEpmTUo8o/TFvrsov+vbbb/K9/tlnn0efG9NyVP3tu5CqqvIKTrsRmbdbR5V+3ByNj4YGVvVuERlXoDRk6P3o06fgQMjP37+0D03GXhagLtuWEJk7/VylFQc0PeCCPW2rereIjCdQ8vf3Vycyr/5ukrPfthYTuYnMXX65ShxVIlNU6kBJ6+DB/3Ds2FEkJibC28sb7Tt0UD3gyHTEJmfjVHiG2m4caA93x6JrZBGR6eOoEpmDUidzZ2Vl4eOPPsAnkz7C0SNHEBcXhx07t2PC+Fcwb97c8t1LqlJ7L9ycdmtfh9NuRGQ4qiS0o0pEpqbUI0qbNm7A1athmPnVLIMpuMOHD+HTTz7GHXd0Ru3atctrP8lI6idJoUkiooJGlcb08kaQB3OVyHSUekQpLCwMffvedVueUvPmLdCiZStcDbtSHvtHRpTIbW0FtA7liBIRGY4qScNc7ajS7G2xPDxkUkodKPkHBCA2Nibf6+LjYrnqzUSEx2fh0o0u4S2CHeBoW+qXDBGZqJF3eMDBRtOdYcneeMQkZ1f1LhGVm1J/6nXu1Bl79uzG4sWLcO3aVSQlJarGtl9/NQO2tnYIDg6uFsUnJddq3bo1mDFjOn7++UfVw45u2n1jtZvgtBsR5cfdyRrD2nuobekH+ft2jiqR6bDIk262pTDrm6+wadPGYt3WWItPpqQk471335bDgLvuultdtn79Orw28XV4e/sU6zFktZ+bmxsSEhLg6uoKU/PGomtY+V+C2v5tdE20qcWpNyq5zIx0RERchr9/CGzt7HkITVBEQhb6fX4W2TmAi70lNk6sC2d7rpCl6q/CCk5Wh+KTv8/9DSmpqZg2bQYcHBzUZd179ERu7o0S1GZOYmht/SQZVm8WrDlGRES38nezwaBWbli6LwFJ6blYsDsOT3X35oGias9sC05mZ2djy5bNuP/+EbogSdjYsIeZ1sXoTEQmanINJInb1lqTg0BElJ8nunlh2f4EyDzFb//G4pFOnrCzYV4jmXnByeoqIiIcGRkZCAmpiUWLFqhVfD7ePujdpw8CA4MKzWmSk1Zq6s0aQyZdFqAuywIQUeFq+djhziYuWH80CTHJOVh+IAEjOmhyl4iqK5MKlCTgOX/+fKG3GTduAmxtbZGenq7O//DDd+jcpQvatW2HQ4cO4pWXX8JHH3+C+vUb5Hv/pUsWq+fR0g+aTM3u83r93VhokoiK4akeXipQEr9sicHQtu6wtuJoNFVfJhUotWjRssj2KVZWmuRCJydn9bN169YYNeoxtd21W3dERkZi+bKleG3iG/nef8jQYRg46F6DZO6AgACYmtzcPOy5UT/J1cESDQKYgEtERWsS5IBOdZ2w42wKwuKysPZIIga0dOOho2rLpAKlBg0aFvu2vr6+cHR0hL+/YZAjeVcyDVcQyWHSz2OSXCdTdDI8HQlpmrIO7Ws7wcqS3wiJqPijShIoiZ+3xKB/C1dYSEdtomrIbLPsZGSpe/ee2L17l276LCkpSU2/NWzUCOZOW41bdOC0GxGVQPvajmgerBmFPh2Rga2nknn8qNoy20BJjHz4YTU69Pxzz+DDD97Diy88i+DgEAwf/gDMnbYsgGChSSIqCRk9errHzdIAP27Ov4sDUXVgUlNvJeXg4IgPPvwYFy9cQFx8HPz8/Apd8WYuMrPzsP+iZkTJ19UatXzY4JKISqZHQ2fU8bXFuchM/HcpDfsvpLJgLVVLZj2ipBVaqxZatWrNIOmGI2FpSMvUFGzvUNuRuQVEVGKWlhZ4sruX7vyPm6N5FKlaYqBEt9mtVz+pA+snEVEp3dPCDYHumsUv206n4MQ1TVkWouqEgRIVWj+pQx0WmiSi0rGxssBjXT1152UFHFF1w0CJDKRm5uLgZU2gFOJlo/s2SERUGkPausPTSVO/bt2RRFyKzuSBpGqFgRIZ+O9iqur+LTiaRERl5WBriUc6a0aVcvOAX7dxVImqFwZKZGCXXv0klgUgovLwQEcPONlpPm6W709AZKLptn4i08NAiQzs1qufJEXjiIjKytXBStccNysnD7/9G8uDStUGAyXSkZYlx2+sSqnvbwdPZ7Mus0VE5WhUF0/YWmvamCzcHa9rkURk7Bgokc6Bi6nI05RPYn4SEZUrHxdr3NfGTbdoZP7OOB5hqhYYKJHOkStpuu1WNR14ZIioXD3ezQva/tq/74hFWmYujzAZPQZKpHPkys1icM2DGSgRUfkK9rTF3c1d1XZcSg6W7I3nISajx0CJlNzcPBwN04woebtYwd+N+UlEVP7025pIqQDpLUlkzBgokXI5JhOJ6Zph8GY1HNjfjYgqRIMAe3Rv4Ky2IxKy8dehBB5pMmoMlEg5EsZpNyKqHE/1uDmq9PPWGDWiTWSsGCiRclgvkbsZ85OIqAK1DnVEm1BNHuT5yEz8fSKJx5uMFgMlum3FW5Mgex4VIqpQT/fw1m3/uDkGedraJERGhoESITM7FyfDM9SRqOVjq6roEhFVpC71ndAgwE5tHw1Lx2699klExoSBEqkgSdoKaBO5iYgqmoWFBZ7SWwH34+ZoHnQySgyUyGDarVkwp92IqHL0beqKYE8bXUPuEzdaKBEZEwZKpIa9tTiiRESVxdrKQlXr1pq9LYYHn4wOAyXSrXizsbJQNU6IiCrLoFZucHfU5EWuPZyIiIQsHnwyKmZbfjk7Oxt/rl6V73UNGjREw0aNYA6kg/fF6Ey13SjQTtfdm4ioMjjYWuKBjh747u9oZOdqesBNuNuPB5+MhtmOKMlS1Pj4eIPTuXNn8dtvsxFxPQLm4tiNtiWC025EVBUe7Oih+5K2eE88UjJy+Icgo2G2I0o2NjZ49LHHDS77/fff4OTkhDvu6ARzrMjNQpNEVBW8XawxsKUrluxLQFJ6LpbsTcCoLp78Y5BRMNsRpVvl5ORg8z9/o1u3HrCz09T2MLeK3M1ZkZuIqsioLjeTuuduj0X2jZIlRFWNgdINBw7sR1xcHPrc2bfQA5aVlYXU1FSDU3WeftSWBnC1t0SIl2aZLhFRZavrZ6eKUIpr8VnYeIxtTcg4mNTU286dOxB5/Xqht+k/YCCsrW//tTdt2oh69eohNDS00PsvXbIYixYtMAicqqvwhGzEJGtyAZoGO6gCcEREVeWxrl7493SKrlRAv2YufF+iKmdSgVJKcrJKyi5Mfv2E5D4H9u/D6NFjinyOIUOHYeCge3XnExMTERAQgOpeaJLTbkRU1TrWcUR9fzucjshQ+ZMHLqWhTahjVe8WmTmTCpSKmjYriOQmSXJ35y5di7yt3E5O+mUGTKIidw3WTyKiqiWj2o919cSb/wtX5+dsi2GgRFWOOUoy7fb3RnTp2g0ODubV50x/xVtT9ngjIiNwd3M3+LpqvsP/fSIZl27UeSOqKmYfKJ04cRzXrl7FnX1KNxpVXcmKEm0NpUB3G7U8l4ioqkk9pYfu8FDbkinx2/bYqt4lMnNmHyjJqrURDzyIuvXqwZyci8xAWpYmX6s5G+ESkREZ3sEDDraaxSXL98cjPqX6pjhQ9Wf2gVKbNm0xfPgDMDdH9Ctys34SERkRNwcrDGnjrrbTs/KwcE/hi3SIKpLZB0rm6sgVvYrczE8iIiPzSBdPWN6oWPLHzlhkSiM4oirAQMlMaStyW1kCjYK44o1Mz9WrYdi7Z3dV7waVUrCnLXo3cVHb0Uk5+PNgIo8lVQkGSmYoNTMXZ69n6KrhOtryZUCm58Tx41iyZHFV7waVwWN6/d7m/Bubbx08oorGT0gzdOJqOnJvvN9w2o2q0rVr17CnHEZ9yutxyLi0rOmIliGasi1nrmdg+xlN1W6iysRAyQzpN8JtxhVvVIVOnjyBxf9baDSPQ8bn0a56o0rbWCqAKh+L55gh/RVvbF1CVUVaB504fgyJiUnYuGG9uqx9h45ISkpUtc2aNG2GM6dPq+r3Xt7eSEiIR4sWLXX3P3ToIFxdXeHh4Znv42ilpaXi5MmTyMjIQJMmTeHiosl7oeqhd2MXBHva4EpsFnacTcHpiHTU92deJVUejiiZ8Yo3qVNSx9euqneHzFRaWhoiIiKQkZGO06dPqVN6errKLfrppx/w1puvY8OGdbh8+RIOHNiHNWv+Mri/nN+/b1+BjyPi4mLxzttvqqbXK5Yvw9iXnkdcXFwV/cZUGlaWFniks2GuElFl4oiSmYlOysa1+Cy13STIQb0JEVUFaSbds1dvrF3zJ557/kWD62JjY/HGm+8gNDRUnV+6dHGpHichIQHvvf8RAgMD1fmJr03Atm1bMUivsTUZv8Ft3PH1higkpudi9cEEvNzXBz6uN3tuElUkjiiZdaFJDl+TcQoMCtIFSWVRo0YNXZAkQkJCEBMTXebHpcrlZGepqnWL7Bypq8RRQao8DJTMzBH9RG4WmiQj5WBv2KDaArePfObmFl2A0Nra5rbu9HnaJZ9UrYzs5AFrK832gt1xqswJUWVgoGRmjoTpVeRm6xKqYra2NsiWIYIiSAJ2rN5IUFZWFi5dvFjix6Hqy9fVBvc0d1PbiWm5qgccUWVgoGRGcnPzcPTGiJKXsxUC3JiiRlWrZs1QVUFbcpBkxVpiYv7Vl1u1ao0rV67g++9mYfXqVfj4ow+QmZlZ4sch0ykV8Nu/scjh6CBVAgZKZuRyTKZKhtSWBZBpCKKqFBwcgtdff0slXWtXqwXVqIH27TsY3E7KA0yZMhVOzs5ITEzAw4+MwvARI1CrVu0SPU6jxk3QoEGDSv0dqfw0DLBHxzqOalvKBfxzIomHlyqcRR5rwpeJfHN1c3NTb9BS08WYrfovAa8vuqa2X+rrg2d6elf1LpGZyMxIR0TEZfj7h8DWjosIqPS2nUrGmNlX1Harmg74fUzZk/6JCsMRJXOtyF2DH1ZEVP10qe+EOr62avu/S2k4fPnm+xpRRWCgZKYr3ppyxRsRVUOSMvBoFy/d+dn/xlTp/pDpY6BkJjKzc3EyPENt1/KxhavDjXW2RETVzICWrmpBithwNAlhsTcT+4nKGwMlM3EqPANZOZr6MayfRETVmZ2NJR7sqClAKQvf5m5nWxOqOAyUzDA/qSnzk4iomnugowfsbTQrd5fuS0BiGutoUcVgoGQmjuoVmpTSAERE1ZmHkzXuba0pQClVupfvT6jqXSITZfaBUk5ODq5du4bz584hKSnJ5Hu82VhZoEGAXVXvDhFRmT10x80ClH/sZAFKqhhmXZr5+PFj+GrmlypYcnV1U5V9u3TthjFjnoOVlekkOyek5eBClCbZsWGAHWytzT4+JiITUNfPDp3qOmHH2RRVgFJqLPVo5FLVu0UmxqwDpa+/molGjRrjhRfHwtLSEpcvXcKECa+gadNm6N69B0zFsRujSYLTbkRkas1yJVAS83bEMVCicmfWQwvSCqFe/foqSBIhNWvCyclJXW5K2AiXiExVtwbOCPa0UdsSMJ29rimDQlRezHpE6aGRj2DlyhVwdHSCu5s7duzcDncPD/To0bPA+0jXcjlppaRovskYcxPO/acjkZOh2c9Qt0yj3lcy3RYmSUnJcHRMhK0da95Q+bqvuTWmr4tT279suoTXB/jzEFOxubi4FNr71KR6vV2+fBkpKcmF3qZBg4a6EaTIyOuYNetrXA0LUzlK0dHRGDnyYfTtd1eB91+4YD4WLVqgO5+amoq16zaW429BRERElaWoXq0mFSjNmf2r6hxemHff+wB2dnZqVOilF59Dixat8PToZ1Ty9pUrl/HG66/hsceeQJ87+xZrRCk3N1cFWB4eHoVGpNWFBH6jn34CP/z4CxwdNV26iceKryv+HxorvmfxWFX0iJJJTb09+tjjxb7t5cuXEBkZqUaPtCvcgoNDVCL33r17CgyUbGxs1Emfs7MzTIW1tbX6/SS6ZqDEY8XXFf8PjR3fs3isKprZJnPLVJuIiYk2uDwmJgaubprriIiIyLyZ1IhSSfj4+KBd+w746cfvVaFJD3cP7Ny5HZcuXcQzzzxb1btHRERERsBsAyUxfvyr2LBhPQ7s34fUtFT4+wXg86nTULNmKMyVTLsNH/7AbdOLxGPF1xX/D40R37N4rCqaSSVzExEREZUns81RIiIiIioKAyUiIiKiAjBQIiIiIiqAWSdzm3MV0pjoaHj7+BRajVSfrAy8HhEBLy8veHh6wlxIixr5vaWgaFG/d1paKi5evHjb5aGhteDg4ABzoa2QX69efVXjpjBSvPXKlSuwtbVFjRo1YG4iIiIQFxeLkBBNn8mCnDlzGtnZ2QaXeXv7qNW7pkxWIUtBSX3ynhUUVPRr5dq1q0hPS0dwSIhZLE4JDw9HfLymjYuWvb09atWqXa73MUcMlMyIfIj/+stP6g3Ezc0dV6+GoU2btnjxpZdVtfKCSMuWpUsWw8/fXwUNXbp0xbPPvaAr1GmKoqOjMPvXX3Ds2FH1gRQefg21a9fBK+MmqKApP1cuX8Hbb72BBg0awMLi5mDtc8+/UKw3dlNw7txZdQwyMzPx08+zCzxW4vDhQ5g+7QvYO9gjLTUVXl7eeOPNt9TxNgexsTF4682JiI+PxwcffqyK3RZk8qeTVAFYbf030aNnL9xZQGFcU/HD99+p/0X910Tz5i0w4oEHC7xPXFwcpkyehGvXrsHZ2UV92XnppZfRpm1bmLJlSxerYsmBgUG6ywICAvDCi2PL9T7miIGSGZHedg8/Mkp90xcyqjRhwjgsX7a0wDee/fv2YfH/FuH9Dz5C48ZN1JvP6xMnIKRmKAYNuhemKioqCr1734kJr05U5+Vb7Vtvvo7f5vyKsS+PK/S+77z7gVmNIOmPqEngc/c9/bFi+bJCbysjTl9M/UxVxh858hE1svTBB+/iq5kzVNBg6qT10Ywvp6tgR/7/imPY/SPQvXsPmJtevfoUGhjd6ttvv5EF3fjxp1/VF0B5/5o27XPMmvU93NzdYcqaNWuOceNfrfD7mBvmKJmR9u076IIk4eXtrb49yDfbgmz6eyOaNGmqgiQRGBiIzl264u9Npt0IuFGjxmjVurXuvHybr127NmIKOVb6w9kXLpxHeno6zMn3332Ldu3bo1nT5kXedveuXer43HffUHVepkbuGzwUR48eUQG9qVu8eBHs7O3Rr5AG3LdKTEzA2bNnkBAfD3MiNe7OnjmjvtgVRUbnpC7evYPv042SD7zxhW77ju0wdZlZWWpUV6Z0JRivqPuYG44omRn5Rzh16qSaGjl6RPOhNGbMcwXeXj7wO3XqbHBZvbr1sHHDejUKYOpz/5IbIh/o58+dw4ED+/Hyy+OLvI8M+9va2SHy+nXcdfc9GDXqMZOephQSOIeFXcHzL7ykXldFkddVQECgQT/BuvXqaa47fx6+vn4wVSdOHMeG9esxdeo0ZGRmFPt+ixYuUMdFpoHlC48ca19fX5i69evW6t6rfHx88eJLYwvMobl48QKkNGCdOnV0l0nAFBwcrF5zpu6/A/tx/XoE4uPiVd7fmGefQ6tWrcv9PuaGgZKZkeDm97m/qQ9/ecOVD/KAwMACb5+clKw6K+uT8xJwydy/u4kPZS9Z/D814nb16lV07HiH7sM8P84uzmqKUoayxelTp9R0kpurG4YMHQZTJbluc+fOwYcfTSp24JycfPvrSttcOik5GaZKFkV8OX0ann3ueTUNVNzRs4dGPoIePXqqgFsWY3z6yceYPm0qPvl0SqFdz6s7mZqV/zsJduQ9a+aM6Zgy+VNM/3JmvtPbyclJ6udtry0XVyQnaa4zVa3btMXIh0fBzc0NOTk5mPvbHEz9fIo6VgV98SjNfcwRp97MjLzhTPpkMr6Y9iVmfvUNdu7Yjjmzfy3w9tbWVsjKzDK4TEajNNeZfpz9+htv4bPPp+H7H35W+Vny4VQQSYjUBkmifoMGKgfl33+3wZR98/VMtG3bTgU/Mlpy5cpldfnZM6cLDATkAz8rS/M60g/iTf11JTluMt0tH/JyrM6dPasuv3zpku645ad37z66UUn5UHvwoZE4ffoUIiMjYcokJ0s7hSarsR57/AlERUWqLyH5sbLSvHYyb3vPyjDp15WQgFJeG0JeK4+MelRt79u7t1zvY45M+5VDhZKVJJJvJMFSQWT58a05TDGxseqNXjsCYA7kd+3d5058/90s9c2ruFNp7m7uheaAmQJ3dw8VRMpIpZCRRrFkyf/QvUdP3H13/9vuI1Mohw4dMrhMe5xMecm7TDVqR3X1g0PpORkTE6P7oCqKdiRXjpmfn/l885fVuqKgXEHta0eOizYAUOdjYlGnTl2YE3mPklV/xcmrLMt9zAEDJTMiQ9fyrUxfRHi4wTC1JIpeC7+G+vUbqH+a5i1aYtvWrQbBwd49u9USXbM7VhHhqtaN9jjISjip86Ktk5TffQ4dPqhq5Jiy1ya+YXD+vwMH8PHHH2Di62/pygPI60dGQGTUTT7AWrRoqcpOSO5X7Rv5JHv27FaBhP6CA1Pz+BNPGZyXEbdnx4zGk089bVAeQI6Vh4en+uDP93V18KB6HZpy2YmMjAw1lWtpeXPi49Chg+pnSEiIQdkTGxtrdSwkd0nqLMmSd20ek+TOSZqBvOZMlaRCSNCtX+YlLCwMMTHRBu8/8oUmIyNdHZvi3ocYKJkVWY4ttYDq1a+vEh5l6b98OOl/0O0/sF9Npfw2dx6cnJwxYMAg/PP33/jii8/Rs2cv9SEoK28+/fQzmLI5s3+BtY2NWvEnCY7Hjx/DqpUr8OhjT+huI0GS1Aya8tlU1K1bD3N/mw1YWKjpN8kb2bJls5paee+9D2HuZJRJjtVLY19R0ykNGzVS03XTp0/FAw+ORFJiIhYumI+HRj6sjre5+2TSRyo/56GHHsaxo0exevVKdO3aTRU9PXnyBFauWI77hgwtdsHY6kiCyK+/mqmmHX39/NT0pIxSdunaTf2/aX337TcqqV2WuEvw+OBDD+OXn3+Eo6MTvDw9sXDhAvU/2bJlK5gqCXhee208+vS5EzVqBKvyJpJfKcdJfzGOHL+LFy6o1Ivi3ocYKJkVqQm0fv1abNy4QX2bCAwIxJczvjIoNiZD+g0bNoKlpZXu/KeTP1O1Xv5cvUpV5p70yRSE1qoFU/bU089g8+Z/VH6RfKP38/XFpEmTDZK5ZfRDjpU2qfSJJ5/W3UfyuKTS9FdfzzKbAopaTs5O6rjo54TIB5hcpj8dMn7Ca1i9aiU2bdyggqPnnn9RBQPmxMbGVh0X+VDXJyO6Mj0ppFCiHNN//t6EqOgo+Hj74M233jHIhzNFwcEheOHFl7Bu7Vrs2r1TjbDJCt07bvkQlxFd/cKmffv2UwHkls3/qP/dTp07Y9CgwTBlMir07rsfYM2aP9Voo5OzM4YOHabSBfTTBOQ936KE9yHAIk+GFoiIiIjoNlz1RkRERFQABkpEREREBWCgRERERFQABkpEREREBWCgRERERFQABkpEREREBWCgRERERFQABkpEVKUVu3fv2qm2paTbv/9uVW10qtr169dV1XpphWEqpH3FyRMnqno3iKodBkpEVGXm/zEPx44dU9tSLX76tC9wJexKse8v95HgKjExsdz2SZpEjx83VlXClorFpkIqMU+ePAkx0dFVvStE1QoDJSKqEpGRkdi4cT3uu2+IOi/98Tp37mLQ5qQo2dnZKri6ejWs3PZr06aN6Nfvbkx8/U089fRomAppP9Sh4x343/8WVfWuEFUrDJSIqMykq/u5c2d15yMiwlXT5cKsW7dG9SuTRq/aQKl9hw5wdXHVBUEyWpScnKyCKnk8/ecQu3fvUj+PHD6sbitNm7Wkz5fs14ED+xETE2NwP/3HliBr166dampKLpOf0THRals605fl8bRkFGfvnt1q9CwtLe22aT65TkbHLlw4j3379qoGpfmJjY1R04HSGDcnJ+e26wt7HiENibdu3YzU1NR8H5+IbnezayURUSnt2b1LBTGTp0xV511d3fDWrNfxxhtvGzQS1ieBT69evXXntVNvH3z4Mdzc3ZGRka7Ot2nbDuHh1xAQEIiTJ46jbdt2eGnsK+o+/x3Yr34eP3EMYWFX4OXtjVatW+O//w5g5ozpquGzg6MjTp08ifuGDMGQIcPU7bWP3bp1G1wLv4ZaobXQ586+2LN7N1JSknHl8mXk5eZCOmGGhoaW+vECAwPx+++/Yf26tajfoCFSU1JUYDR+wqto3LiJuu+xY0cxZ/YvCA4JAfKgmglLIPTii2PRuUtX3fGZO3cO/vpzteruLsdKAqW33n4HLjcCy6KeRzRo0FD9PHLkMDp06MjXO1ExMFAiojJr3aYt1q1bi7jYWDVC5OjoiDvu6IxVq1filVfG33b7rKwsFdjUCA4u8rFtrK3x5ZdfqY7mFy9cwPjxL+PewfehZs1QjHn2eWzZshkjRjyIRo0aq9tLMvgXUz/Dy6+MV0GVCAsLw2uvjkPzZi0MAjdJIJ8x42sVnIiWLVup23W8o5MuCCrL423bthVbt2zGjJnf6Drcr169CjNnfIlvZn2n69IuI1G9e9+Jnj17qfOLFi7AvHm/6wKlzf/8jdWrVuLjSZ+iXr366jIJTDMzMgGX4j+P/JRg7+zZMwyUiIqJgRIRlZlModnY2GDvvr3o27efukw+sE+ePJ7v7ZOSElVQ4ezkXORj39m3n+6DPrRWLTg5OeHa1asqUMrPzp07VKCSlZmJHTu2q8vkuSSAO3rsiEFg06/fXbqgpiBleTxJCA+pGapGwvJu3M/WxgZRUZFqxEdGnIS9vT169Oipu1+Tpk2xaNECFVDKcd28+R906txFFySJOnXqlvh5hLOzM5ISk4o87kSkwUCJiMrM1tZWBUuSH6MNlI4ePQJfP/98b29v76B+pmekF/nY8sGuTwKHzKysAm8fGRWpC3D01aldB54emnwoLY9bzpf340lulRybW+8rSet5ebm6846OTipHS/c7WtuoYEdyn+T3jYqOQpMmTQvex2I+jzbXysFBc/yJqGgMlIio3KbfJNdGPoglx+bw4UN46+13872tTM25u7sj8vr1cj/6jg6OalRn3PhXi76xRcU+nqOjA+rWrY/Rz4xBWTg5OiEpqeBRoJI8j4wy9e5zZ5n2h8iccNUbEZWLNm3aIjMzE39v2ohvvp6JO+7opJKbCyIjUBJQlYWMtkgQI8+r1bJVK8TFxanRLX0ZGRkqF6ikyvJ4LVu1xo4d/94W5Ny6aq4oLVq2VKNF8pxa8junpaWW6HnCw8MRHx+P5s1blOj5icwZR5SIqFz4+PggJKQmfv75R5Wk/fwLLxZ6e1kV9tmUT9UHvkwblYZMV9WqVVslOkvStawAk1VvgwcPwbRpU3HXXfcgqEYNXI+IwK5dOzB+wsTbpvKKIqvMSvt4UiPq4H8HMPG1Cejbr58anZJE6osXL+Czz6cVex/k+aVswOsTX0Xv3n3UqjdJ4B4/4TU4ODgW+3m2b9+mAlQ/P78SHQMic8ZAiYjKzcCBg9SU28OPjFIf4IVp2rQZQmrWVCu6+va767aCk9bWNuq8s7OLwf3at+8AX19f3XkJFtavX6uW8Lt7eKhA6ZFRj6qfMgokCc5BQTXw4YeTdDWbCnps0aJlKwTXMFyNV9rHk2Mw6ZMpKqg5deokrCyt0LBRYzw9+uYUmQQt7du3N7ifi6urejxtYrgksE+e/LlK2pYASMonjH15HPz9/Yv9PJIYvmH9erzw4kuF/l2IyJBFnmQMEhFVgUuXLmLzP//g0cce5/GvYBJASTD5wAMP8VgTlQADJSIiIqICMJmbiIiIqAAMlIiIiIgKwECJiIiIqAAMlIiIiIgKwECJiIiIqAAMlIiIiIgKwECJiIiIqAAMlIiIiIgKwECJiIiIqAAMlIiIiIgKwECJiIiICPn7P/PxK6Pg8U8VAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def profile(nu):\n", " res = minimize(lambda t: -loglik_fast(mz, None, nu=nu, p=1 / (1 + np.exp(-t[0])), kernel=kernel),\n", " x0=[fit.x[1]], method=\"Nelder-Mead\", options=dict(xatol=1e-6, fatol=1e-8))\n", " return -res.fun\n", "\n", "\n", "grid_nu = np.linspace(2.8, 5.8, 31)\n", "prof = np.array([profile(v) for v in grid_nu]) - ll_hat\n", "cut = -chi2.ppf(0.95, 1) / 2\n", "inside = grid_nu[prof >= cut]\n", "print(f\"95% profile interval for nu: about [{inside.min():.2f}, {inside.max():.2f}] (grid step 0.1)\")\n", "\n", "fig, ax = plt.subplots(figsize=(6.4, 3.6))\n", "ax.plot(grid_nu, prof, color=accent, lw=2)\n", "ax.axhline(cut, color=ink2, lw=1, ls=\"--\")\n", "ax.text(grid_nu[0], cut + 0.15, \"95% cut-off\", color=ink2, fontsize=9)\n", "ax.axvline(nu_true, color=grid, lw=1)\n", "ax.text(nu_true, -7.6, \" truth\", color=ink2, fontsize=9)\n", "ax.set_xlabel(\"ν (interference)\", color=ink2)\n", "ax.set_ylabel(\"profile log-likelihood\\nminus maximum\", color=ink2)\n", "ax.set_title(\"Profile likelihood for ν\", color=ink, loc=\"left\")\n", "ax.set_ylim(-8, 0.6)\n", "for side in (\"top\", \"right\"):\n", " ax.spines[side].set_visible(False)\n", "ax.tick_params(colors=ink2)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "5c9d87b8", "metadata": {}, "source": [ "## Notes\n", "\n", "- **Per-meiosis parameters.** `nu` and `p` may also be arrays with one value per meiosis, for\n", " example `nu = np.exp(b0 + b1 * age)`. The same recipe then fits regression coefficients. With\n", " many distinct `nu` values the Chebyshev kernel falls back to the exact series.\n", "- **Phase-unknown data.** For pooled crossovers from two meioses, replace `loglik_fast` with\n", " `loglik_quartet`; nothing else changes.\n", "- **Coordinates.** Use genetic-map positions in Morgans. Physical positions would confound\n", " recombination hotspots with interference." ] } ], "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.12.11" } }, "nbformat": 4, "nbformat_minor": 5 }