{ "cells": [ { "cell_type": "markdown", "id": "a6d69e48-7399-4589-a8d6-8954689fa723", "metadata": {}, "source": [ "# Visualization of phylodynamics inference results" ] }, { "cell_type": "markdown", "id": "17b4aa73-4560-4b19-9899-039d0e4cf1a9", "metadata": {}, "source": [ "In this notebook, we will introduce how to visualize scPhyloX inference and analysis results" ] }, { "cell_type": "code", "execution_count": 1, "id": "04a9bf42-2713-44d3-8f35-62522befe79c", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] } ], "source": [ "import scPhyloX as spx\n", "import numpy as np\n", "import pandas as pd\n", "import gzip\n", "import arviz as az\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from matplotlib.gridspec import GridSpec\n", "import matplotlib.lines as mlines\n", "import pickle\n", "from scipy.integrate import quad, solve_ivp\n", "from scipy.stats import poisson, chisquare\n", "import os" ] }, { "cell_type": "code", "execution_count": 2, "id": "af891e8b-8cb4-466c-bf35-af357745c42a", "metadata": {}, "outputs": [], "source": [ "plt.rcParams['pdf.fonttype'] = 42\n", "plt.rcParams['font.size'] = 12" ] }, { "cell_type": "code", "execution_count": 3, "id": "b11c552b-bfbf-4e67-98e7-880034061aa1", "metadata": {}, "outputs": [], "source": [ "os.chdir('../..')" ] }, { "cell_type": "markdown", "id": "630207d5-3454-4aa1-9b6e-6c839867bc9a", "metadata": {}, "source": [ "First, we import the data and scPhyloX inference results. The point estimates of each parameter will be saved in `theta_h`" ] }, { "cell_type": "code", "execution_count": 4, "id": "9e721594-2215-4347-905f-2773c317985b", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/wangkun/miniconda3/lib/python3.9/site-packages/arviz/utils.py:184: NumbaDeprecationWarning: \u001b[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\u001b[0m\n", " numba_fn = numba.jit(**self.kwargs)(self.function)\n" ] } ], "source": [ "idata = pickle.load(open('datasets/simulation/overshoot/para_inf.pkl', 'rb'))\n", "idata_bl = pickle.load(open('datasets/simulation/overshoot/mutrate_inf.pkl', 'rb'))\n", "theta_h = az.summary(idata).loc['ax,bx,r,d,k,t0'.split(',')]['mean'].to_numpy()" ] }, { "cell_type": "code", "execution_count": 5, "id": "80e99ed3-4b5a-4a49-9fc7-bde735df67b4", "metadata": {}, "outputs": [], "source": [ "with gzip.open('datasets/simulation/overshoot/character_matrix.csv.gz', 'rb') as f:\n", " charater_matrix = pd.read_csv(f, index_col=0)\n", "charater_matrix = charater_matrix.to_numpy()\n", "with gzip.open('datasets/simulation/overshoot/simulation_data.csv.gz', 'rb') as f:\n", " ground_truth = pd.read_csv(f, index_col=0)\n", "ground_truth = ground_truth.to_numpy()\n", "ground_truth = ground_truth[np.arange(0, ground_truth.shape[0], 1000)]\n", "time = ground_truth[:, 0]\n", "cell_number = ground_truth[:, 1:]" ] }, { "cell_type": "code", "execution_count": 6, "id": "694e599d-1db5-48b8-be84-8744136fcdcf", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 125250/125250 [00:04<00:00, 25575.91it/s]\n" ] } ], "source": [ "mutnum = spx.data_factory.get_mutnum(charater_matrix)\n", "branch_len = spx.data_factory.get_branchlen(charater_matrix)" ] }, { "cell_type": "markdown", "id": "cd696cca-63f8-4e1f-85d6-fc614382ee7a", "metadata": {}, "source": [ "Here, we can inference the number of stem and non-stem cells changes over time in different generations, that is, the tissue development history could be revealed." ] }, { "cell_type": "code", "execution_count": 7, "id": "4532992c-06cd-4f6e-b190-a1b6278a3675", "metadata": {}, "outputs": [], "source": [ "T = 35\n", "n_stem = np.array([[spx.est_tissue.ncyc(i, j, 100, *theta_h) for j in range(T)] for i in range(100)])\n", "n_nonstem = np.array([[spx.est_tissue.nnc(i, j, 100, *theta_h) for j in range(T)] for i in range(100)])" ] }, { "cell_type": "code", "execution_count": 8, "id": "eb8eec57-d1d2-4467-b1cc-29b3ad3432d2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Cell number')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAE4CAYAAAAjCop8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABqbElEQVR4nO2dd3hUVfrHP/feqcmkF1pCQq+hKdKkKCJFd1UQXUVWRNG1txUXxACCyqq4i13kJ7C6FnTFgrJKEXVFuooIFnoLSUgv0+/5/THJwJBJJSHtfJ5nnjDnnHvve+8M3znlPe+rCCEEEolEIkGtbwMkEomkoSAFUSKRSEqQgiiRSCQlSEGUSCSSEqQgSiQSSQlSECUSiaQEKYgSiURSghREiUQiKUEKokQikZQgBTEIBQUFTJ8+nUsvvZS4uDgURWHOnDm1dv5Zs2ahKAo9e/astXNKJJKzRwpiELKysli8eDFOp5Mrr7yyVs/9ww8/8Mwzz9CiRYtaPa9EIjl7DPVtQEMkKSmJnJwcFEXh5MmTLFmypFbO6/F4uOmmm7jtttv48ccfOXnyZK2cVyKR1A6yhxgERVFQFKVKbd99910GDRpEaGgoNpuN0aNH8/333wdtu2DBArKzs3n88cdr01yJRFJLSEE8C5544gmuu+46unfvzooVK3jjjTcoKChg6NCh7N69O6Dt7t27mT9/Pi+//DI2m62eLJZIJBUhh8w15MiRI8yePZu77rqL5557zl8+atQoOnXqxNy5c3n33XcB0HWdqVOnMn78eMaNG1dfJkskkkqQglhDPv/8czweD3/+85/xeDz+covFwvDhw/nyyy/9Zc8++yy///47H3/8cX2YKpFIqogUxBqSnp4OQP/+/YPWq6pvNuLw4cOkpqayYMECTCYTubm5gG+BRdd1cnNzMZvNWK3Wc2K3RCIpHymINSQ2NhaA999/n6SkpHLb7d+/H7vdzr333su9995bpj4qKop7772Xf/7zn3VlqkQiqSJSEGvI6NGjMRgM7Nu3jwkTJpTbrk+fPgHD51Luu+8+8vLyWLp0KQkJCXVpqkQiqSJSEMth9erVFBUVUVBQAPhWid9//30Axo0bR3JyMo899hiPPPII+/fvZ8yYMURFRZGens6WLVsIDQ1l7ty5REZGMmLEiDLnj4yMxOPxBK2TSCT1gyKTTAUnOTmZQ4cOBa07cOAAycnJAHz00UcsWrSI7du343Q6admyJf379+cvf/kLI0eOLPf8I0aM4OTJk+zatasuzJdIJDVACqJEIpGUIB2zJRKJpAQpiBKJRFKCFMQShBDk5+cjZxAkkuaLFMQSCgoKiIiI8K8qSySS5ocURIlEIilBCqJEIpGUIAVRIpFISpCCKJFIJCVIQZRIJJIS5F7mauD1enG73fVthqQW0TQNg8FQ5ZQRzRGhexFuO3jsCHcxwm1HeJygexC6B3QvCA/C6wHdA5S6rilQ+lwVxfdeNaCohlN/NQOoRhTVgGIwg2ZGMZhRDBYwmFGUc9tnk4JYRQoLCzl69Kj0U2yChISE0KpVK0wmU32bUq8IIRCuIoQjB92ei3AWIDzFPvGrzetUo62imXxCabSiGEJQjCEoRmvJKwRFM9aqbVIQq4DX6+Xo0aOEhIT48zRLGj9CCFwuF5mZmRw4cIBOnTr5A/s2B4TuRdiz0e056PYchCPX1+NrQAivC+F1gTO4f7CiGvwiiTEEQ0xHFK3mP2xSEKuA2+1GCEFcXJyMbN3EsFqtGI1GDh06hMvlwmKx1LdJdYrQvehFGegFaeiFGQjhrW+TzgqhexDOfHDmA2CI7XJW56v3n8P169czdepUunbtSmhoKG3atOGKK65g+/btVTo+IyODKVOmEBsbS0hICIMGDWLdunV1YqvsGTZNmnqvUOhevAVpuI/vwLV3De7jO/AWpDV6MTwTRTOhqNpZnaPevwkvv/wyBw8e5N577+Wzzz5j0aJFZGRkMHDgQNavX1/hsU6nk5EjR7Ju3ToWLVrERx99RIsWLRgzZgxfffXVOboDiaRhIlxFuNN34drXdEXwdBRjyNmfo77jIWZkZBAfHx9QVlhYSMeOHenZsydr164t99iXXnqJO++8k40bNzJo0CDAl7ypd+/e2Gw2Nm/eXGU78vPziYiIIC8vj/Dw8IA6h8PBgQMHaNeuXZMfUjVHmtrnqzvz8Wbtw1uQRvWWMBo3WlgrjK37ndU56r2HeKYYAthsNrp3786RI0cqPHblypV06dLFL4YABoOBG264gS1btnDs2LFat1ciaajojlzcx7bhOvgN3oLjNCcxhNrpITbIRZW8vDx27NjBxRdfXGG7Xbt2MXTo0DLlvXr1AuDnn3+mTZs2dWJjYyEjI4NHH32U1atXk56eTlRUFL1792bOnDkMGjQIRVFYuXIlV155ZX2bKqkhuiMXT+av6MUnz8HVFBSjxbeyqxpBNYCq+X0LUQ0lvoM+MfYNQEXJW93nu1jqr6i7favaXjfC6wavEyH0mltmPPsFzwYpiHfeeSdFRUU88sgjFbbLysoiOjq6THlpWVZWVrnHOp1OnM5T/lX5+fk1tLZhM2HCBNxuN8uXL6d9+/akp6ezbt06srOz69s0yVkidC+ek7/izTlQ6+dWVAOKNQrVEoFiDC0RwdBadZYWHhfCmYfuyEM4ctEdeaiWKAyxncHrQHicvleJMzilTuHlzIM2yR7io48+yr///W+ef/55zjvvvErbV7TyW1Hdk08+ydy5c2tkY2MhNzeX//3vf2zYsIHhw4cDkJSUxAUXXADgT5R11VVX+esOHjwIwCeffMKcOXP4+eefad26NTfeeCOPPPIIBoPvK6MoCq+88gqffPIJ69evJykpiddff524uDhuueUWtm7dSq9evXjzzTfp0KHDub3xZoBedBJ3+k6fUNQCqjkMxRKFao1EsUShmELPyqtCCB3hLPQLnXDkojvzEfZcdEcuwpkX1HZDfHdMbfoBtvLP7XH5HMbdp+2ccRc3PUGcO3cu8+fP5/HHH+euu+6qtH1MTEzQXmBp7ydY77GUGTNm8MADD/jf5+fnk5iYWAOrGy42mw2bzcaHH37IwIEDMZvNAfVbt24lPj6epUuXMmbMGDTN57Lw+eefc8MNN/Dcc88xdOhQ9u3bx6233grA7Nmz/cfPmzePZ599lmeffZaHH36Y66+/nvbt2zNjxgzatm3L1KlTueuuu1i9evW5u+kmjvC68WTuxpt39KzPpVoiUMNaoYW1qhUxKcW+6z086bugBsNf3ZFXaRvFYEIxmMASWQPrKqbeF1VKmTt3LnPmzGHOnDnMnDmzSsekpKTw008/lSkvLevZs2e5x5rNZsLDwwNeTQ2DwcCyZctYvnw5kZGRDBkyhJkzZ7Jz504A4uLiAF+O6JYtW/rfP/744/ztb3/jxhtvpH379owaNYp58+bx6quvBpz/pptu4pprrqFz5848/PDDHDx4kEmTJjF69Gi6devGvffey4YNG87pPTdlvAUncB386qzEUDWHY4jrirndCExJF2KI7lBtMRS6F+EqKrdeUbQaiSGAcOTW6LjaokEI4rx585gzZw6zZs0K6IFUxlVXXcUvv/wS4F7j8Xh48803GTBgAK1bt64LcxsVEyZM4Pjx43z88ceMHj2aDRs20K9fP5YtW1buMdu3b+exxx7z9zBtNhvTpk0jLS2N4uJif7vSxSuAFi1aAL4fqdPLHA5Hk52fPVcIoeM+sRP38e013les2VpgajsYU/JQnwiaQqtxfYFefBLXkc3Yf/w3RV89iePXT8ttr1Sz56YYQ1DDWmGI64qhRa96jRdQ70PmhQsXkpqaypgxY7jsssvYtGlTQP3AgQMBuPnmm1m+fDn79u0jKSkJgKlTp/Liiy8yceJEFixYQHx8PC+99BK//vprhf6LzQ2LxcKoUaMYNWoUqamp3HLLLcyePZspU6YEba/rOnPnzmX8+PFBz1WK0XhqY33pfFOwMl2v+cphc0d4XbiPbUe312QRTEELa4UW0wHVXL0RkHAX48nejzdrL57sfWV6bt7sfQghgs4zqpaI00zQfAszlghUS2TJ3wgUS6S//Gz2Htc29S6In3zyCQD//e9/+e9//1umvvTXwuv14vV6A349zGYz69atY/r06dx9990UFxfTp08fVq9e7V9EkJSle/fufPjhh4BPwLzewFW7fv368euvv9KxY8d6sE5Siu4qwn10K8Jd/vA0KIqCFp6IIbp9tXqCAJ7s/bj2f4k392CF7YS7GL3wBFpYqzJ1WmwXQvrf5hM7k61RbXmtd0Gs6hzTsmXLgg7zWrRowfLly2vXqCZCVlYWEydOZOrUqfTq1YuwsDC2bdvGU089xRVXXAH4VprXrVvHkCFDMJvNREVFkZqayuWXX05iYiITJ05EVVV27tzJTz/9xPz58+v5rpoHenE27mPbEHr14m+qlkgMLXuhmsOqdZy34ATOvV/gzfq96sfkHw0qiKrZBubyV4kbMvUuiJK6w2azMWDAAP7xj3+wb98+3G43iYmJTJs2zb9wtXDhQh544AFee+012rRpw8GDBxk9ejSrVq3iscce46mnnsJoNNK1a1duueWWer6j5oE37yju9J1Qjbk0RdHQ4rqgRSZXq0emO/Jw7VuHO+37yq9hsqFFd8AQ0xEtun21h+GNgXrfy9xQkHuZmy8N6fP1nPwVT9beah2jhsRibJlS7dVi3Z5D0XfPlUS5DnZiA1pUOwzRHdCiO6DaWjSq4W9NkD3Es6DYXvNgmiajisEQfJHfbvciguxDDbHKj6sp407fhTf3UJXbK6oBQ3wPtIiEGl1PtUahRbUrM0xWTDZM7S/G2KpPrUekbujI/2FnwWv/2l/jY0cMiad3z8igdW+sOIjdUXZ70r23da7x9SQNG0/mr9UTQ2MopoT+1V40ORNzx0spLhVEzYgpaSimtkN8js/NECmIEkk948k5gCe76sNk1RqNsc15VXZX8eYdQQ1PCDrc1cJaYmx9HqgapvYXo56lwDZ2pCBKJPWIN/8YnozdVW6vhSdgaJlSpQALQvfg/P0L3Ee+w5JyLcYWwXdumbtd0eTnBqtKg9ipIjn3TJky5ZyE/JozZw59+vSp8+uUh6Iofp/Lhoa3MAP3iR+r3N4Q2wVjq95VEkNvYQbFW17FfeQ7AFx71yD0cqLESDH0I3uIZ8G0P7ev8bEmY/lf6snXJAddVKlNFi1a1GBTqjaHGI26PQfP8R1Vcq1RFBVDy95o4ZVvRRVC4D6+DeevnwWsHuv2bNxHt2BqO6iCoyVSEM+Culr1tVrPLlFOVYiIiKi8kaRO0J0Fvh0oVchvoigqxoT+qCGxlbYV7mIcez4KPgRXtAa1Ra6hIofMTZz333+flJQUrFYrMTExXHLJJRQVFZUZMo8YMYK7776b++67j6ioKFq0aMHixYspKiripptuIiwsjA4dOgSE8lq2bBmRkZEB1/vwww8rHIJt3bqVUaNGERsbS0REBMOHD2fHjh3++tNjNCqK4n8Pvm2e5513HhaLhfbt2zN37lw8nlO9oN9//51hw4ZhsVjo3r07a9asqdlDq0OE24776JYq7kBRMLTqWyUx9OQcoGjTi0HFUA2NI2TA7RjbVB5ftLkjBbEJk5aWxnXXXcfUqVPZs2cPGzZsYPz48eUOlZcvX05sbCxbtmzh7rvv5vbbb2fixIkMHjyYHTt2MHr0aCZPnhwQ8aa6FBQUcOONN/LNN9+wadMmOnXqxLhx4ygo8CUi37p1KwBLly4lLS3N/740RuM999zD7t27efXVV1m2bBmPP/444AsgMX78eDRNY9OmTbzyyis8/PDDNbazLhBCx338e4THUaX2xpYpaGEtKz6n7sW5by327a/78hOfeY6ECwi54HY0W4sa2dzckEPmJkxaWhoej4fx48f7IwSdHp7rTHr37s2sWbMAXwDdBQsWEBsby7Rp0wBITU3l5ZdfZufOnf4oRNXlzDw5r776KlFRUXz11VdcfvnlZWI0lnJ6jEaA9u3bM2/ePKZPn87s2bNZu3Yte/bs4eDBgyQk+ByVn3jiCcaOHVsjO+sC78nf0R05VWpriO2CFlFxwGLdnoNj13t488omY1OMVizdr8IQ161GtjZXpCA2YXr37s3IkSNJSUlh9OjRXHrppVx99dVERUUFbX96fENN04iJiSkT3xB8iatqSkZGBqmpqaxfv5709HS8Xi/FxcUcPny4wuO2b9/O1q1b/T1C8EVAcjgcFBcXs2fPHtq2besXQyAgG2N9oxefrLKvoSGqHYaYiiMNeTJ/wbHrfYS3bHxELaodlh5Xo1qa3l7jukYKYhNG0zTWrFnDxo0b+eKLL3j++ed55JFHys1XfXosQ/Ct9lYU31BV1TLDb7e74rmxKVOmkJmZyT//+U+SkpIwm80MGjQIl8tV4XGVxWgMNg3QUNxJhMeFO+2HKrXVwtugVdKr82Ttxb7zHThzUUZRMXcYiTHpwlpLBNXckILYxFEUhSFDhjBkyBBSU1NJSkpi5cqVtXLuuLg4CgoKKCoqIjTUt8Phhx9+qPCYb775hpdeeolx48YBcOTIEU6eDEyfWZMYjd27d+fw4cMcP37cHyn9u+++q8lt1SpCCNwnfqxSpGs1NB5Dy16VCrk77fsyYqhao7H0nFjjfc0SH1IQmzCbN29m3bp1XHrppcTHx7N582YyMzPp1q2bP6/K2TBgwABCQkKYOXMmd999N1u2bKkwNQFAx44deeONNzj//PPJz8/noYcewmoNzKdbkxiNl1xyCV26dOHPf/4zCxcuJD8/v9I0tucCb85B9KLKpxgUkw1j635V6tlZuo/HqWj+kF2GFilYul2BYjBXcqSkMmS/ugkTHh7O119/zbhx4+jcuTOzZs1i4cKFtbbQEB0dzZtvvslnn31GSkoKb7/9NnPmzKnwmNdff52cnBz69u3L5MmTueeee4iPjw9os3DhQtasWUNiYiJ9+/YF8MdoXLNmDf3792fgwIE8++yz/sUiVVVZuXIlTqeTCy64gFtuuSVgvrE+0B15eE7uqbSdoqgYW/VFUavmf6qoGubuV2FsO8gnhj0nSjGsJWQ8xBJkPMTmS118vkL34Dr4DcJduYuSIb4Hhqjk6l9DCEDI+cJaRA6ZJZI6wJP+c5XEULO1qJEYQumiUcNYOGoqyJ8WiaSW8Rak4c2vPHeyYrBgaNGr3HrhceHY8xF6BTmQJbWLFESJpBYRuhdPRuXzhqD4IlKXE4hV6F7sO9/CfWwb9m2voTvyatdQSVCkIEoktYg363eEx15pO0NMB9SQmKB1Qui+HSjZ+wDQi7Mo3roYvehk0PaS2kMKokRSS+iuIjw5laeVUK1RaDGdgtYJIXD+thpPxs9nnNyNEHptmCmpACmIEkkt4cn4udL4hopq8LnYlLMy7D70P9xHNgUWakasff6MZosPeoyk9pCCKJHUAt6CE+hFmZW2M8R1RTFag9a5037AufeLwEJFw9p7ktyBco6QgiiRnCW+hZTK86KolgjUiLZB6zxZe3HsLrul0tJjPIboDmdto6RqSEGUSM4Sb/beqi2kxPcIuk/ZW5CGfedbcMYcobnTGIwty3fLkdQ+UhAlkrNAuIrwZle+kKJFJKJay4Zd0+052L//F3gDowQZ2w7ClDSk1uyUVA0piE2cKVOmoCgKCxYsCCivLNR/XXPw4EEURak0Ok5Dx52xu9LVX0U1YojtWqZcuIuxf78c4SoMKDfE98DcqeEEtm1OSEFsBlgsFv7+97+Tk1O1aM2SquEtzKhSJBtDXJcyDtjC68L+w5voxVkB5VpUMpaeVzeYWI7NDSmIzYBLLrmEli1b8uSTT5bb5j//+Q89evTAbDaTnJzMwoULA+qTk5N54oknmDp1KmFhYbRt25bFixdXeN2cnBwmTZpEXFwcVquVTp06sXTpUgDatWsHQN++fVEUhREjRviPW7p0Kd26dcNisdC1a1deeuklf11pz3LFihUMHToUq9VK//79+e2339i6dSvnn38+NpuNMWPGkJlZ+apvTfEtpPxcaTvVHB50IcX5+xdlQv+rofFYe12PosoQA/WFFMRmgKZpPPHEEzz//PMcPVp2j+327du55ppr+NOf/sRPP/3EnDlzePTRR8vENly4cCHnn38+33//PXfccQe33347v/zyS7nXffTRR9m9ezerV69mz549vPzyy8TG+jLIbdmyBYC1a9eSlpbGBx98AMBrr73GI488wuOPP86ePXt44oknePTRR1m+fHnAuWfPns2sWbPYsWMHBoOB6667junTp7No0SK++eYb9u3bR2pq6tk8tgrR849WLZJNi55Be3um9hehhrfxv1fM4Vj7/rlclxzJOUJIhBBC5OXlCUDk5eWVqbPb7WL37t3CbrfXg2Vnx4033iiuuOIKIYQQAwcOFFOnThVCCLFy5UpR+vFff/31YtSoUQHHPfTQQ6J79+7+90lJSeKGG27wv9d1XcTHx4uXX3653Gv/4Q9/EDfddFPQugMHDghAfP/99wHliYmJ4q233goomzdvnhg0aFDAcUuWLPHXv/322wIQ69at85c9+eSTokuXLuXadjrV/Xx1r0c49q4V9l9WVfhyHf+h4vO4naJox3JR8OV84SlIr9K1JXWL7CE2I/7+97+zfPlydu8O9Jnbs2cPQ4YErmgOGTKE33//PSCU/+lJqBRFoWXLlv6EU2PHjsVms2Gz2ejRowcAt99+O++88w59+vRh+vTpbNy4sUL7MjMzOXLkCDfffLP/XDabjfnz57Nv376AtqfbUpr86syEWGeTDKsi9PyjlaYSVVQDhriyCykBbQwmrL0nYT1/mtyF0kCQkxXNiGHDhjF69GhmzpzJlClT/OVCiDLDOhFkC1qwJFSlCaeWLFmC3W4PaDd27FgOHTrEp59+ytq1axk5ciR33nknzzzzTFD7Ss/12muvMWDAgIA6TQuMJh0s+dWZZaXnq02E7sWTVXn2PC22c5WiWCuqJsWwASEFsZmxYMEC+vTpQ+fOnf1l3bt353//+19Au40bN9K5c+cyQlQebdq0CVoeFxfHlClTmDJlCkOHDuWhhx7imWeewWTyrbqe3gNt0aIFbdq0Yf/+/UyaNKm6t3ZOqFLv0GBFK2dHiqRhIwWxmZGSksKkSZN4/vnn/WUPPvgg/fv3Z968eVx77bV89913vPDCCwGruzUhNTWV8847jx49euB0Olm1ahXduvlSbMbHx2O1Wvnvf/9LQkICFouFiIgI5syZwz333EN4eDhjx47F6XSybds2cnJyeOCBB87KnrOlqr1DQ0yHgPwoQgi8J39Fi+0i3WkaOHIOsRkyb968gCFxv379WLFiBe+88w49e/YkNTWVxx57LGBYXRNMJhMzZsygV69eDBs2DE3TeOeddwAwGAw899xzvPrqq7Ru3ZorrrgCgFtuuYUlS5awbNkyUlJSGD58OMuWLfO76dQn3rwjVeodquGBgRg86Tux//hv7Nv/T8Y0bODUe5KpgoIC5s2bxw8//MD333/PyZMnmT17dqXZ2wCWLVvGTTfdFLQuLS2Nli1bVtkOmWSq+VKVz1foXlwHNlQqiMYWPdEik/zvdVcRxd89d8pFR9Ewdx6LKXFAOWeQ1Cf1PmTOyspi8eLF9O7dmyuvvJIlS5ZU+xxLly6la9fAFb2YmODRiGuTs8l1oWgmFM0YtE64igj2K6WaQmt8PcnZUeXeYURiQJnz188C/RWFF8UcVhcmSmqBehfEpKQkcnJyUBSFkydP1kgQe/bsyfnnn18H1lVM0dcLKm9UDuYul5fbSyg6vUdxGmGXzKvx9SQ1R+hevNlVnDs8LfCrJ/NXPOk7A9vEd8cY373WbZTUDvUuiHKSWdLQ8fUOnRW2ObN3KDwuHL98fEYbC+Yul9eJjZLaoVqLKna7neuvv76Mi0Z9c/nll6NpGtHR0YwfP55du3bVt0mSJkJNe4euAxsQzvyANubOY1HlcLlBUy1BtFqtfPTRR3Xi8FoTWrZsySOPPMKSJUv48ssvmTdvHlu3bmXgwIH8+OOPFR7rdDrJz88PeDUnpkyZwpVXXlnn15kzZw59+vSp8+uUh6IofPjhhzU+via9Q73oJK7DgbtytKj2GFr1rbEdknNDtYfMffr0YdeuXQwbNqwu7KkWY8aMYcyYMf73w4YN47LLLiMlJYXU1FQ++uijco998sknmTt37lldP3TY32p8rKIFz8cLEDronqCLKrXJokWLgu5GaQgoisLKlSvPiWBXhBB6lYK/nt47FELg+HUViFMO5ygalq5/kNNDjYBqC+KCBQuYPHkyPXr0YPjw4XVh01mRnJzMhRdeyKZNmypsN2PGjABH3/z8fBITEys4oix1teqrmEKp6/86ERERdXyFxo9emF5paoAze4eezD3+fMqlmNoOQg2NrRMbJbVLtR2z77jjDgoLC7n44ouJjY0lJSWFXr16+V+9e/euCzurhRACVa341sxmM+Hh4QGvpsj7779PSkoKVquVmJgYLrnkEoqKisoMmUeMGMHdd9/NfffdR1RUFC1atGDx4sUUFRVx0003ERYWRocOHVi9erX/mGXLlhEZGRlwvcoicW/dupVRo0YRGxtLREQEw4cPZ8eOHf765ORkAK666ioURfG/B/jkk08477zzsFgstG/fnrlz5+LxePz1v//+O8OGDcNisdC9e3fWrFlTs4dWgjfnQKVtDDEdT/UOvS6cv30WUK+YwjC1u+is7JCcO6otiDExMaSkpDBs2DBSUlKIjY0lJibG/4qOjq4LO6vMgQMH+Pbbbxk4cGC92tEQSEtL47rrrmPq1Kns2bOHDRs2MH78+HKHysuXLyc2NpYtW7Zw9913c/vttzNx4kQGDx7Mjh07GD16NJMnT6a4uPI4gOVRUFDAjTfeyDfffMOmTZvo1KkT48aNo6CgAPAJJvh8S9PS0vzvP//8c2644Qbuuecedu/ezauvvsqyZct4/PHHAV9giPHjx6NpGps2beKVV17h4YcfrrGduiMP3V5xhHFf7/DUrhTXwW8QjryANubOY8pEy5Y0XKo9ZN6wYUOtG7F69WqKior8/yl2797N+++/D8C4ceMICQnh5ptvZvny5ezbt4+kJN9OgEsuuYRhw4bRq1cvwsPD+emnn3jqqadQFIV586TPXlpaGh6Ph/Hjx/uf2ekhss6kd+/ezJo1C/BNKSxYsIDY2FimTZsG+PYmv/zyy+zcubPGPzgXX3xxwPtXX32VqKgovvrqKy6//HLi4uIAiIyMDNhp9Pjjj/O3v/2NG2+8EYD27dszb948pk+fzuzZs1m7di179uzh4MGDJCT4ROqJJ55g7Nia5SapSu9Qi0r29w714ixcB78pU29oUf7zljQ86t0PEXxx8w4dOuR//9577/Hee+8Bvh5fcnIyXq8Xr9cb0LtJSUnh3Xff5ZlnnsFutxMfH8/FF1/Mo48+GhDNpbnSu3dvRo4cSUpKCqNHj+bSSy/l6quvJiqqbPY3CIwxqGmafzRQSmncwbOJM5iRkUFqairr168nPT0dr9dLcXExhw8frvC47du3s3XrVn+PEHyRchwOB8XFxezZs4e2bdv6xRBg0KBBNbJReBx4C45X2EZRtIDk8a7D3wUupKBg7nK5XEhpZNRIEDMzM3n22WfZsGEDJ0+e5MMPP6RHjx68+uqrXHDBBfTtWz33goMHD1baZtmyZWVC2v/jH/+o1nWaG5qmsWbNGjZu3MgXX3zB888/zyOPPMLmzZuDtg8W7zBY3MFStytVVcsMv93uwHSaZzJlyhQyMzP55z//SVJSEmazmUGDBuFyuSo8Ttd15s6dy/jx48vUWSyWoNMANRUjb+5hqGQFXo1ICPAUMHcei2qNxrV/HcLrwth2IJqtRY2uL6k/qi2IBw4cYMiQIeTl5dG7d2/279+P0+nz09q5cyebNm3yJxKS1D+KojBkyBCGDBlCamoqSUlJrFy5slbOHRcXR0FBAUVFRYSG+lbcK0sr+s033/DSSy8xbtw4AI4cOcLJk4ERYIxGY0CcRPBF5Pn111/p2LFj0PN2796dw4cPc/z4cVq3bg3Ad999V+17EroXb+6hSttpkckB7xVVw5Q0GEPLFFwHNmCWCymNkmovqkyfPp3IyEh+//13vv7664Bf5gsvvJBvv/22Vg2U1JzNmzfzxBNPsG3bNg4fPswHH3xAZmamPybh2TJgwABCQkKYOXMme/fu5a233irTiz+Tjh078sYbb7Bnzx42b97MpEmTsFoDEyslJyezbt06Tpw44U+dmpqayr/+9S/mzJnDzz//zJ49e3j33Xf9c56XXHIJXbp04c9//jM//vgj33zzDY888ki170kvSEN4K+6tqiGxqGZb8DpzmM/nUCaLapRUWxDXrVvH7Nmzad26dZkhSatWrTh+vOK5F8m5Izw8nK+//ppx48bRuXNnZs2axcKFC2u80HAm0dHRvPnmm3z22WekpKTw9ttvVxq27fXXXycnJ4e+ffsyefJk7rnnHuLjA0PoL1y4kDVr1pCYmOiffhk9ejSrVq1izZo19O/fn4EDB/Lss8/6F4tUVWXlypU4nU4uuOACbrnlloD5xqpStcWU+o/NKKkbqh0PMSQkhI8++ohRo0bh9XoxGo1s27aNfv368emnn3Ldddc1ym1wMh5i86X0801qGYWasb3CtooxFFO74XKxpIlS7R5ily5dWLt2bdC6r7/+mp49e561URJJfeAtOFZpG5+rjYLQvejOgnNgleRcUu1FlWnTpvHAAw/QunVrfyIgl8vF+++/z0svvcQLL7xQ60ZKJHWOEL7w/hV0ERTV4He18ZzYieOXjzG26Y8p+UJUc9Pc6dTcqFEKgVtvvZUlS5agqiq6rvvdL6ZNm8Yrr7xSF3bWOXLI3HxxOBzs3/sbbbTjmFVvue20qHYY47sjhE7xxufQ7Vm+CkXD3PFSTEmDz5HFkrqiRn6IixcvZurUqaxatYqMjAxiY2O5/PLLGTxYfiEkjQ8hBOgeqCTjqqEkV4on/adTYgggvKghdZ+yQlL31HinysCBA+V+YUnTwFuxMzmAZmuBYgpFCB3Xga8C6tSwVmixcmdUU6DGgvjdd9/x5ZdfkpWVRUxMDCNGjJA9REmjROhVEMSoZAA8GbvRizID6sztL5Krzk2Eagui3W7nT3/6E6tWrQpwylYUhXHjxrFixYoyjrYSSUNF6N4z9iCXRTWHoYbEIoTAtf/LwDpbS7TYruUcKWls1GinyurVq5k/fz4HDhzAbrdz4MAB5s2bx+eff8706dPrwk6JpE6obFcKnNqm58ncg14UGNjC1G6E7B02IardQ3z33Xd59NFHmTFjhr8sKSmJmTNn4na7eeGFF3j++edr1UiJpC4QQkAlgqioGmp4GwDcZ4T3UkPjMMiUok2KavcQi4uLy50rHDJkCHZ7xSHXJZIGg+6GSrLXqGGtUVQNb+4hvPlHA+rkjpWmR7UFceDAgf4oxmeydetWLrjggrM2SlJ7ZGRkcNttt9G2bVvMZjMtW7Zk9OjR/kgwZ5uVrjFTpeFySb4U16HAoCWKJQJDvNyV1dSo9pD5ueee47LLLiMsLIzrr7+eqKgocnJy+Pe//83ixYtZtWpVXdgpqSETJkzA7XazfPly2rdvT3p6OuvWrSM7O7u+TatXhO71+R5WgGoOQ7VGoRdn4cncE1BnShyMolbiuChpdFRpp0pYWFjA0MDlcvkDgRoMBn+iH6PRiNlsJi8vL+h5GjJNcadKbm4uUVFRbNiwIWiGxOTk5IBI5UlJSf5gvZ988ok/1Fbr1q258cYbeeSRRzAYfL+hiqLwyiuv8Mknn7B+/XqSkpJ4/fXXiYuL45ZbbmHr1q306tWLN998kw4dOpyT+60Owu1AeB0AOJwuDh46ShtjRsBOFUNcdwzR7XD88gnuo1v85YpmJnToX1EMjee7IKkaVeohTpgwQc6VNEJsNhs2m40PP/yQgQMHYjabA+q3bt1KfHw8S5cuZcyYMWiar8dTmtDpueeeY+jQoezbt49bb70VgNmzZ/uPnzdvHs8++yzPPvssDz/8MNdffz3t27dnxowZtG3blqlTp3LXXXcFZOprCAghKh8uKwpaeBuEuxj38R0BVcaE86UYNlWERAghRF5engBEXl5emTq73S52794t7HZ7PVh2drz//vsiKipKWCwWMXjwYDFjxgzx448/+usBsXLlyoBjhg4dKp544omAsjfeeEO0atUq4LhZs2b533/33XcCEP/3f//nL3v77beFxWKp5Ts6e3SPS3jtOf5XUW66+PnH7SJ392ph/2WVsP+ySriObRdCCOE6tl3kr5l16rU2VXjtufV8B5K6otqLKpLGxYQJEzh+/Dgff/wxo0ePZsOGDfTr16/CyNbbt2/nscce8/cwbTYb06ZNIy0tLSAF6elJqUoTUJ2ZlMrhcDS4+JjVWUwxtu5HSP/bSrLnKRhapKBaIurYQkl9UaOte06nk7Vr13Lo0CEcDkdAnaIo3H///bVinKR2sFgsjBo1ilGjRpGamsott9zC7NmzmTJlStD2lSV0KiVYAqqKklI1BITQK11MUQxWlJBY/3stIgFryjXoHUfVtXmSeqbagrh9+3b+8Ic/kJ6eXm6mMymIDZvu3bv7XW1qktCpUeOt3PdQi0gIOmeuWoOnb5U0HaotiHfccQfh4eG88sordOvWDZPJVPlBknohKyuLiRMnMnXqVHr16kVYWBjbtm3jqaee4oorrgBOJXQaMmQIZrOZqKgoUlNTufzyy0lMTGTixImoqsrOnTv56aefmD9/fj3f1dlRneGypPlRbUH8+eefeeutt/jjH/9YF/Y0WIQAp1MHpeJAAOcSk1FFVctf/bfZbAwYMIB//OMf7Nu3D7fbTWJiItOmTWPmzJmAL6HTAw88wGuvvUabNm04ePCgP6HTY489xlNPPYXRaKRr167ccsst5+rW6gSheyoP5BASKzPmNWOqHTG7e/fuLFiwoMkJYmV+iL/9dgAPcRgNNY6YVut06hCGxSydg6uK7i4Ounf5dD9EW0Jv1JBYUA3S1awZUqNoN88884w/Ob1E0hgQQlQaCFZRjaih8Tj3fkHx5pdwH9uOqELwWEnTodrdnSlTpnDw4EE6dOjAiBEjiI6ODqhXFIVFixbVmoESSa1QhcUUJTQedA+e498jvE4cez5E2fs5lh4TMMR2OTd2SuqVagvip59+ypNPPonb7eatt94qUy8FUdIQEXoVFlPCW+FO+wHhPTX6EW47amh8XZomaUBUe8j80EMP0a9fP3744QecTie6rge8znThkEjqG5/vYSXfS0VBMYbiPro5oNgQ21W62zQjqt1DPHjwICtXrgzYpSBpPkyZMoXc3NxGFTLsphtvJCcnm5Xv/bv8RqoBb95RxBn5UoyJA+rYOklDoto9xK5duza4rViS2ufgwYMoisIPP/wQUL5o0aIKt/3VFlOmTOHKK6+slXMJUdlOGQVF0fCc+DGgVA2JQYtueJF6JHVHtQVx3rx5zJ8/nxMnTtSFPZIGTkREBJGRkfVtRpURupfKFlNQNAQCb/a+gGJj4kDpetPMqPaQ+dVXXyUnJ4eOHTvSp0+foKvMH330Ua0Z2FAwGCCpTQhmi7nyxucIk7Hy3zMhBE8//TSvvPIKaWlpdO7cmUcffZSrr76anJwc7rrrLr744gsKCwtJSEhg5syZ3HTTTbRr1w6Avn37AjB8+HA2bNhQZsg8YsQIUlJS0DSN5cuXYzKZmDdvHpMmTeKuu+7i/fffJz4+nhdeeIGxY8cC4PV6ufXWW1m/fj0nTpygbdu23HHHHdx7770AzJkzh+XLlwOn9kN/+eWXjBgxgmPHjvHAAw/wxRdfoKoqF154IYsWLSI5Odl/7oceeojXX38dTdOYOuXGoFtMT0fRDOBxnlFmwtiqbxU+BUlTotqCuHPnTjRNIy4ujmPHjnHs2LGA+qb6i6ooYDarjc4RetasWXzwwQe8/PLLdOrUia+//pobbriBuLg43nvvPXbv3s3q1auJjY1l7969/pw4W7Zs4YILLmDt2rX06NGjwi2ay5cvZ/r06WzZsoV3332X22+/nQ8//JCrrrqKmTNn8o9//IPJkydz+PBhQkJC0HWdhIQEVqxYQWxsLBs3buTWW2+lVatWXHPNNfz1r39lz5495Ofns3TpUgCio6MpLi7moosuYujQoXz99dcYDAbmz5/PmDFj2LlzJyaTiYULF/L666/zf//3f3Tr1o2FTy/gw48/5aLhQ8uxXgHFUGZLn6FVXxRDw/nxk5wj6jP2WEOiKcZDLCwsFBaLRWzcuDGg/OabbxbXXXed+MMf/iBuuummoMceOHBAAOL7778PKL/xxhvFFVdc4X8/fPhwceGFF/rfezweERoaKiZPnuwvS0tLE4D47rvvyrX1jjvuEBMmTCj3OkII8X//93+iS5cuQtd1f5nT6RRWq1V8/vnnQgghWrVqJRYsWCCE8MU9dBZkioQ2rcUfLx8XEAPR/3IViaK8bLFrx0Zxct0Cf9xDb2FGubZKmi4NZx+apNbZvXs3DoeDUaMCw1a5XC769u3LnDlzmDBhAjt27ODSSy/lyiuvLDejYkWc7nGgaRoxMTFl4iKCL+FVKa+88gpLlizh0KFD2O12XC4Xffr0qfA627dvZ+/evYSFhQWUOxwO9u3bR15eHmlpaQwaNAgAobsxGAyc169vucNmRTUhPLkBZVpUO9TQuErvW9L0qLYgHj58uNI2bdu2rZExktqlNA7hp59+Sps2bQLqzGYziYmJHDp0iE8//ZS1a9cycuRI7rzzTp555plqXef0GIjgmzapKC7iihUruP/++1m4cCGDBg0iLCyMp59+ms2bA30Ag93Peeedx7//XdZ9Ji4uUMBEFbbq+dYURUk60tPup03/So6TNFWqLYjJycmVzhNK5+yGQffu3TGbzRw+fDhokinwCcmUKVOYMmUKQ4cO5aGHHuKZZ57xzxnWxWf5zTffMHjwYO644w5/2b59gSu8JpMpaJzGd999l/j4+DIBOEpp1aoVmzZtYuiQgYDA4/Gw4/sf6Nund5m2imZEuIsDy4xWDPHdanhnksZOtQXx9ddfLyOIJ0+e5OOPP+bo0aPMmjWrWucrKChg3rx5/PDDD3z//fecPHmS2bNnM2fOnCodn5GRwfTp01m1ahXFxcX07t2b+fPnM3LkyGrZ0RQJCwvjr3/9K/fffz+6rnPhhReSn5/Pxo0bsdls7Nu3j/POO48ePXrgdDpZtWoV3br5xCA+Ph6r1cp///tfEhISsFgsRETUTuj8jh078q9//YvPP/+cdu3a8cYbb7B161b/yjb4fng///xzfv31V2JiYoiIiGDSpEk8/fTTXHHFFTz22GMkJCRw+PBhPvjgAx566CESEhK49957WbBgAR3aJdKtcwf+8dxL5OaV4zerGRGO3IAiQ6s+KKqcSWquVNsPccqUKdx4440BrwcffJCvvvqKfv36ceTIkWqdLysri8WLF+N0OqvtiOt0Ohk5ciTr1q1j0aJFfPTRR7Ro0YIxY8bw1VdfVetcTZV58+aRmprKk08+Sbdu3Rg9ejSffPIJ7dq1w2QyMWPGDHr16sWwYcPQNI133nkH8KWXfe6553j11Vdp3bq1P6BsbfCXv/yF8ePHc+211zJgwACysrICeosA06ZNo0uXLpx//vnExcXx7bffEhISwtdff03btm0ZP3483bp1Y+rUqdjtdn+P8cEHH2Ty5MlMveUvDBlxKWE2G1f+8bKyRigaKCqKKcz37xKMrc+vtfuUNEJqc4Vm1apVIjExsVrH6LruXzXMzMwUgJg9e3aVjn3xxRcFELCK6na7Rffu3cUFF1xQLTua4ipzc0V3O4KvKJ/20t0Of/vi4mLx866dIu/3r+rRaklDoFaz7nk8HnJzc6t1jKIoNfZdXLlyJV26dPGvKoKvZ3PDDTewZcuWMj6SkuZB5TEMFdACF30U1YAp4YK6NUzS4KkVQXS73Wzfvp3Zs2fTu3fZyeu6YteuXUGDTJSW/fzzz+fMFknDQOheEBVn1UPVUBSZgVdSlmrPHquqWm6PLioqis8///ysjaoqWVlZZbYOAv6yrKysco91Op0BUb9lwIomQhUiXCuaTIwmCU61BTE1NbWMIFosFpKTkxk3blwZp9m6pqLhdkV1Tz75JHPnzq0LkyT1SOWBYBVQjZW0kTRXqi2IVXWHORfExMQE7QVmZ2cDBO09ljJjxgweeOAB//v8/HwSE2X6ycaML6teJaG+1BLfQ9WIYpA9RUkgjdrhKiUlhZ9++qlMeWlZz549yz3WbDZjNsvN+02JqiSEUlQNvfikLyCYakQxhSJE0wxIIqk+NRLE//3vf7z11lv+fainoygK69atqxXjKuOqq67ijjvuYPPmzQwY4Its7PF4ePPNNxkwYACtW7c+J3ZI6h9Rxa16AVFtdDfCkQeGyLo0TdKIqLYgLl26lJtvvpno6Gg6d+5cppclqpfmGYDVq1dTVFREQUEB4AtK8P777wMwbtw4QkJCuPnmm1m+fDn79u0jKSkJgKlTp/Liiy8yceJEFixYQHx8PC+99BK//vora9eurbYdkkaM7gEqGS5rRoQzL6BIMVrlirPET7UF8amnnuKaa65h+fLltTbkvP322zl06JD//Xvvvcd7770HwIEDB0hOTsbr9eL1egME12w2s27dOqZPn87dd99NcXExffr0YfXq1eXu3ZU0TaqUP1mIMsmmFGMoeCpLMSBpLiiiml26kJAQPv74Yy655JK6sqleyM/PJyIigry8vDKBAxwOBwcOHKBdu3ZYLJZ6slBSHkIIhDOfClMFKBp4XQiP41SZakQNjcPpdMrPVwLUwDG7W7dupKen14UtkjpgypQpKIrCggULAso//PDDeo1uXl4SqxqhV56EHkULFENAMYU02QjvkppRbUF84oknWLBggdwW14iwWCz8/e9/Jycnp75NqROqNFzWz2yjoBhC6sQeSeOl2oL44osvkpeXR+fOnRkyZAh//OMfA161GRVFUjtccskltGzZkieffLLcNv/5z3/o0aMHZrOZ5ORkFi5cGFCfnJzME088wdSpUwkLC6Nt27YsXry4wuvm5OQwadIk4uLisFqtdOrUyZ8j5fQkVoqiMGLECP9xS5cupVu3blgsFrp27cpLL73kryvtWa5YsYKhQ4ditVoZMGgov/2+l63bdnDBkIsIj01g7B+vJjPzpO8gRUO4z/CGMFpRVLmYIjmD6kaDSEpKEsnJyeW+2rVrV7vhJ84RTTXaTWlukg8++EBYLBZx5MgRIYQQK1euFKUf/7Zt24SqquKxxx4Tv/76q1i6dKmwWq1i6dKl/vMkJSWJ6Oho8eKLL4rff/9dPPnkk0JVVbFnz55yr33nnXeKPn36iK1bt4oDBw6INWvWiI8//lgIIcSWLVsEINauXSvS0tJEVlaWEEKIxYsXi1atWon//Oc/Yv/+/eI///mPiI6OFsuWLRNCnMr10rVrV/Hf//5X/LzzBzHwgv6iX9/eYsSwC8XX61aLbd99JTp2aC9um3aTL7qNI0948o4GvHS3029nY/58JbWLTDJVQlMXRCGEGDhwoJg6daoQIlAQr7/+ejFq1KiA4x566CHRvXt3//ukpCRxww03+N/rui7i4+PFyy+/XO61a5LEKjExUbz11lsBZfPmzRODBg0KOG7JkiVCCCG8jgLx7+VLBCDWrP7IH97ricdSRZfOnYTXnis8hZmBgliQHpCoqjF/vpLaRY4ZmhF///vfWb58Obt37w4o37NnD0OGDAkoGzJkCL///ntAGP/TIwspikLLli39iaPGjh2LzWbDZrPRo0cPwOdO9c4779CnTx+mT5/Oxo0bK7QvMzOTI0eOcPPNN/vPZbPZmD9/fpkUA7169ULoOggPLVrEA5DSo7u/Pr5FPBmZmYAC3jNyLptC5WKKJCiNeuuepHoMGzaM0aNHM3PmTKZMmeIvF0KUEQgRxBsrWDKp0sRRS5Ys8e9aKm03duzYaiWxKj3Xa6+95t95VIqmBebDNhqNUBLIodT2MxNb6boefDHFaA16fYlECmIzY8GCBfTp04fOnTv7y7p3787//ve/gHYbN26kc+fOZYSoPM7M6ldKdZJYtWjRgjZt2rB//34mTZpU6TWrsrosPMF6h3JgJAmOFMRmRkpKCpMmTeL555/3lz344IP079+fefPmce211/Ldd9/xwgsvBKzu1oTU1NRqJ7GaM2cO99xzD+Hh4YwdOxan08m2bdvIyckJiE7kCwRbhYyAiuLboVL61ihdbSTlIwWxqggd4SpEV6rg83aOUIwhKGrVenCnM2/ePFasWOF/369fP1asWEFqairz5s2jVatWPPbYYwHD6ppQmsTq4MGDWK1Whg4dWiaJ1WOPPUZqaipDhw5lw4YN3HLLLYSEhPD0008zffp0QkNDSUlJ4b777gs8eVV8D1FQbS0QbjvCVYiiGlA0GQtRUj7V3rrXVKls697B338mQTmCwdhw/kOZkoehms9tQN6GgBAC4Syg4mAOCoo53D+/KIQAIYL6HsqtmZJSZA9R0vjQvVQlss3pC0WKoviGzxJJBVRJEC+++OIqn/BcxkOUNE8qTxMAikwTIKkBVRJEXder7LclR+CSukRUMRAsqhz8SKpPlb41GzZsqGMzJJIqUllkGwEYDNLxWlIjpEOWpFpMmTKFK6+8st6uX5nvodDdCHsOuj0H4XEihKh3myWNBzmukATl4MGDtGvXju+//54+ffr4yxctWnROpkWmTJlCbm4uH374ob9MCL0kVUAFlAimcBcj3MVyV4qkWlRJECtKTn8miqLg8VTypZU0WiIiIurv4t7Khsu6z1n79J0omnSjkVSdKglisOT0zQ2hWdASLsTUgFKXVmXXhRCCp59+mldeeYW0tDQ6d+7Mo48+ytVXX01OTg533XUXX3zxBYWFhSQkJDBz5kxuuummgHiFAMOHD2fDhg1lem4jRowgJSUFTdNYvnw5JpOJefPmMWnSJO666y7ef/994uPjeeGFFxg7dizg26536623sn79ek6cOEHbtm254447uPfeewFf7u/ly5f77rHke/fll18ybNB5HDt2nAcfnsWadetRVZUhgwfyz2cWkJzUFuF14/XqTE+dx7I330VTVabefItc6JNUmSoJYkNKTl9vKCqKyYZqblw9jlmzZvHBBx/w8ssv06lTJ77++mtuuOEG4uLieO+999i9ezerV68mNjaWvXv3+gM0bNmyhQsuuIC1a9fSo0cP/97jYCxfvpzp06ezZcsW3n33XW6//XY+/PBDrrrqKmbOnMk//vEPJk+ezOHDhwkJCUHXdRISElixYgWxsbFs3LiRW2+9lVatWnHNNdfw17/+lT179pCfn+8PKBsVGUFxUQ4jx/yRC4cMYsOaTzEYDDy+4BnG/fFqftjyP4x4ePaFxSx7811ee/5puvfsxT9eXMLKlSur5TomacbUT9SxhkdTjIdYWFgoLBaL2LhxY0D5zTffLK677roaxSs8Pb6iEEIMHz5cXHjhhf73Ho9HhIaGismTJ/vL0tLSBCC+++67cm294447xIQJE8q9ju4qFq+98rzo0rmT8BRn++Me2vPShdVqFas/ek94cg+LVi1biCfmzCgJAusQbrdbJCQkBJzrTBrr5yupfWq0yvzLL79w3XXX0apVK0wmEzt27ABg7ty5fPnll7Wn1pKzYvfu3TgcDkaNGhUQX/Bf//oX+/btq3a8wvI4PU6ipmnExMSQkpLiL2vRogWAP3YiwCuvvML5559PXFwcNpuN1157jcOHDwc9vxAC4XWz4/sf2LtvPxFxiYTHJhAem0Bs6/Y4HA727dtLXn4haSfSGXjBeaBqoJkwGAycf/75NbovSfOj2qvMP/zwA0OHDiUsLIwRI0YEBAkoLCzklVde4aKLLqpVIyU1ozS+4KefflomPJfZbCYxMbFa8QrLI1icxDNjE55uz4oVK7j//vtZuHAhgwYNIiwsjKeffprNmzeXcyO+JPS6rnNe3z68seyMXC4CYsMD53YVo8yoJ6k+1RbEv/3tb/Tq1Ys1a9ZgMpl49913/XUXXHAB//nPf2rVQEnN6d69O2azmcOHDzN8+PCgbaoTr7C2+Oabbxg8eDB33HGHv+zMiNgmk8l/7VLfw759erPi/ZXEx8UGBOAQXjd4nKAotGrZgs1bdzBitC/ZmcfjYfv27fTr16/W70PS9Ki2IH777be8+eabhISElPnP0qJFC06cOFFrxknOjrCwMP76179y//33o+s6F154Ifn5+WzcuBGbzca+ffuqHa+wNujYsSP/+te/+Pzzz2nXrh1vvPEGW7du9a9sgy/L3+eff84vv/xCtM1EREQYk/40kYX/eJ6rJk5izqMzSEhow+EjR1n5wUoevOc2Etq05u6/TOWpf75El5Tz6datG88++yy5ubm1Yrek6VPtOUQhRLkrjjk5OZgbkFuKxBf7MDU1lSeffJJu3boxevRoPvnkE9q1a+ePV9irVy+GDRuGpmll4hW++uqrtG7dulbTy/7lL39h/PjxXHvttQwYMICsrKyA3iLAtGnT6NKlC/3796dFYge+/W4zISEhbFjzKYmJCVx93Z/p0WcAt9x2F8X2YsLDfGHQHrj7NiZPup4pU6b4h+NXXXVVrdkuadpUOx7i4MGD6dmzJ4sXL8br9WI0Gtm2bRv9+vXj9ttv5/fff2ft2rV1ZW+dUVk8RBkvr37QXYUV7k4Rboev3j9fqKCGtaxWmgD5+UpKqfaQ+d577+X6668nNDSUyZMnA3D48GHWr1/P66+/zvvvv1/rRkqaJ0KvZKuewBfs4TTx8y2myC36kppRbUG89tpr2bdvH3PmzOG5554DYMKECRgMBubOncsf/vCHWjdS0kypJO6h8JatV0yhdWWNpBlQo+AOM2fO5M9//jOff/456enpxMbGMnr0aJKSkmrbPkkzprLINopqAM2EcBcDAjSzzJkiOStqHO0mISGBm2++uTZtkUj8CN1TeVY9gxnVFIqwhPsi28go2ZKzpEqTLTk5OUyYMIFVq1aV22bVqlVMmDCBrKysWjNO0nwJNhw+E0XzeTsoiopqsqEYpIeD5OyokiAuWbKEH3/8kTFjxpTbZsyYMfz000+8+OKLtWZcQ6OaC/KSGiKEXoU0AUqtpQko3UEjkVTpG/XOO+8wbdo0DIbymxsMBqZNm8a7775LampqrRnYEDAafRncMjMziYuLk1vC6hjhcSG8zoobqSZUKmlT2XWEwOVykZmZiaqqFUb0kTQPqiSIv/32W5U2yPfr14958+adtVENDU3TSEhI4OjRoxw8eLC+zWnyCI8DKuuNG8y15l4TEhJC27ZtUYPkbJY0L6okiB6Pp8wG/mAYjUbc7sqGOo0Tm81Gp06dmuz9NRR0ey7uEz+UW+8tzAAhsHT9A6qp8gC5laFpGgaZlEpSQpUEsVWrVuzevZthw4ZV2O7nn3+mZcuWtWJYQ0TTNDRNq28zmjTu7DRUtZzVZSFwndwJKHi3HkBt1QdT8oWoIbHn1EZJ06VKY4Thw4fz0ksvVdg7crvdvPzyyzL0l6TGCI8Lb2FaufV68UmEqxAMZhBe3Me3oxedPIcWSpo6VRLE+++/n19++YWrrrqK48ePl6k/fvw4V155Jb/++iv3339/rRspaR54849WOHfozd4HBot/eKtaY9Biu5wr8yTNgCoJYq9evXjxxRf94ZoGDx7MpEmTmDRpEoMHD6Zdu3Z88cUXvPjiiwGRkqtKYWEh9913H61bt8ZisdCnTx9/1JWKWLZsGYqiBH3JMGSND29u8IjZAMJZiF6YHpBW1Nh2kJz7k9QqVXbkmjZtGj179uSJJ57gyy+/ZNOmTYBvhW7MmDHMmDGDgQMH1siI8ePHs3XrVhYsWEDnzp156623uO6669B1neuvv77S45cuXUrXrl0DymJiYmpki6R+0ItPItxF5dZ7c/aDZvJt1wMUgwVjq77nyjxJM6Fanq2DBg3ik08+Qdd1Tp70zd3ExsaelbvCZ599xpo1a/wiCHDRRRdx6NAhHnroIa699tpKFzJ69uwp82Y0cirqHeJ14c09jGK2+YuMbfqjGKTfoKR2qZGSqapKfHw88fHxZ+27tXLlSmw2GxMnTgwov+mmmzh+/Hj5eTYkTQbhceAtLH+Kw5t7iNLgDQAoKsbEAefGOEmzot49UXft2kW3bt3K7IIpzeS2a9euSs9x+eWXo2ka0dHRjB8/vkrHOJ1O8vPzA16S+sGbV8FiitDxZO8Hg9U/X2iI74FqqZ10BhLJ6dS7IGZlZREdHV2mvLSsomARLVu25JFHHmHJkiV8+eWXzJs3j61btzJw4EB+/PHHCq/75JNPEhER4X8lJiae3Y1IaoQQAm9e+cNlb94R8DgCFlNMbQefC9MkzZDa2R1/llS0UlhR3ZgxYwICTgwbNozLLruMlJQUUlNT+eijj8o9dsaMGTzwwAP+9/n5+VIU6wG9KBPhtgevFDqek7+VLKb45pG1yCS0iIRzaKGkOVHvghgTExO0F5idnQ0QtPdYEcnJyVx44YX+VfDyMJvNMiFWA8Cbe7D8urwj4C5GOW14bGo3os5tkjRf6n3InJKSwp49e/B4AnNn/PTTT4BvBbm6CCHkRv1GgO7IQy/KDF4pdLwnfwNF8y+maOEJaNEdzqGFkuZGvavGVVddRWFhYZkE98uXL6d169YMGFC91cQDBw7w7bff1tgnUnLu8GbvK78u76gvCrbx1M4UU/uLpCO2pE6p9yHz2LFjGTVqFLfffjv5+fl07NiRt99+m//+97+8+eabfh/Em2++meXLl7Nv3z5/7pZLLrmEYcOG0atXL8LDw/npp5946qmnUBSlSYYha0rozkK8BeXvW1atkajhrfG63Xi8JvTQBLC0w+TyYjKqUhgldUK9CyLABx98wCOPPEJqairZ2dl07dqVt99+mz/96U/+Nl6vF6/XGxC1OiUlhXfffZdnnnkGu91OfHw8F198MY8++iidO3euj1uRVJEze4dCFzhdOh6Pjscj8HjNOMMupdCQgPfkr2imJNTDvp0siqJgNCqYTComo4bJpGK1aISGaFIoJWdFtRPVN1UqSlQvqV2E247zwJcgBG63TlGxh2K7t0yKBnv0YISh6mlFDQaViHAjkRFGQqwN4rde0siQ3xrJOcedtQ97sZuiYi9OV/DYh15zi2qJIYDHo5OV7SQr24nZrBEVYSQywoTRWO9T5ZJGghREyTnD49E5mVmAfe9v6N6KU4y6Q5IB8HoF+YVuCgo9eDw6LeIshIYE/9oePlaMy6kTHm4kItyI0+nlRIYDm81AXIwFW6j8uksqRn5DJHWOEILsHBcnMhxo+b9hrEAMdR1ynBEczIP8/CyKigPb2kKN5QriiXQHhUUeOOZ7bzarRIYbiYgwcSLdQYs4C61aWLFaZdRzSXCkIErqFIfTy7E0O8XFHtDdWOxHy7QRAoqKXHiz9pLvDuNXbxIOHEHPp+vBp7y9XuETw9NwOnXSM52kZ/qy82maQlyMmc4dwujWOQyLRX79JYHIb4SkTtB1QWaWk8yTTv9iidF+BITntDaQm+ciO9uFRT+JTbFjQdBK+ZVM0Y5C4gClzHmDkV9QefIvr1dwIsPBiQwH237IpmM7G+f1jiY8vPIEapLmgRRESa1TVOzhWJodp/O04a7uwVB8BABRIoRZOS68XoGChxB8WzWLRCQmxU4EJygU8QAoCthCDVgtGhZL8OGuxayR3DaE/HwPeQVuvN6KnSeK7V527s5j1y95dGxnY8iAWMLDZHzF5o4UREmtIYSvB3Yyq2wCeYPjGIpwoXvhwOFCPJ5TghVKNqoi8GDCTQgGXBDVmT6tI2kRb6FVvAVriAGhC1xuHZfL93K6dJwuLy6XjtWq0T7pVADZwmIPeXlu8vJd5OZ5cDiDz1vqOvy2r5CQEAPn9Y7CFip7i80ZKYiSWsHt1jl8rNg3V3gmwoux+BAAqgYhVoN/iKviJkTJQ0GgGQx0th0lpn0PIvr1Cnods7lsD7HY7iE3z01evhuPRwfAFmLAFmKgTStf2DCH00t6ppPjJ+zY7YHiGGLViI40ceBQEZERJlrGW6SrTjNFCqLkrCks8nDkWLFfjM7EWHwIRT/Va4yPM1NY5EHXBWFkEm7MJzJUxxJhRFENhHYdVa3rh1gNhFgNtGphoaDQQ06ui4JCT4Cjt8WskZQQQlJCCLl5Lo6fsJOe6UQISG57yt8xN89FQaGbFnEWoqNMcudLM0MKouSsyMxykp5hLzfgteIpxlh0IKBM0xTi48y4ivLpaDyE0aCgmMN9W/LanI8aUr2Qb/5rKQrhYUbCw4y43ToZJ53k5DrL2BYZYSIywkT7ZN8KeIs4S0C91ys4fsJOTp6LI0eLiY+z0K1zOKoqxbGpIwVRUiO8XsHRtGLy84Ov7h44VERMjJk4z6/AqZ6jAoSGGGjVwor32C/ohYovxJfBApqx1uIdGo0qbVpZiY0xkZHpJDfPVaaNxazRIdkW5GgfGZlOdu7OA/L4YVcOI4bE06ZVSK3YJ2mYSEGUVBuH08vhI8VBt915PIKff80jK9tF7omjxLXNBM3XszIaVCIjjJhMGnpxFnpJYinFFIqiKJjaXohqLl+gaoLZpJHYJoTYGDPpGQ4KCit3zyll74FC/79PZrl4/+OjdGpv48KBcYSHycWXpogUREm1yC9wc+RYcVB/wKJiDzt/zsPu8KLgJcKzl2NpHtomhBAeZiQs1ODzoRECT8Zu30GqAQwWFGMopqQhdWa31aKR3DaUoiIPx9PtOBwVbx3ML3CTk1u2V/n7/kL2HSzkvN7RnNcnCrNJ7nppSsilNEmVych0cOhIUVAxzDjpZNv32dhLhCZaOYoRJ3a7F7NJI8xm9IkhoBeeQNh9foeKOczXO2w3HMVQ9ykdQkMNdGxno2ULa4VzguFhRvr0jCQkpKzg6Tps/T6b5W8f5Mefcyv1eZQ0HmT4rxJk+K/y0XXB0ePF5AWZL9QF7D9QyOFjxf4yI8Ukqd9j0CClewTRUacJndBxHdiAcBaAwYJqiUC1RBIy+D5/Iqlzhcutk5ZuL3ceFHzbCo+l2dl/KNB38nQiI4xcOCCO9smhclW6kSMFsQQpiMFxuXQOHS0KOsR0OL38vCePvIJA38M2yi5iQwvomxJVxm/Qm3sIT9oPoCgo1lgUVcXScyLGlsH9Ds8F+QVu0tLtuFzB3YYA3B7BgUOFHD1eToZAoHVLKwPOiyGxjVUKYyNFziFKyqWoyMPhcvwLs3Nc7Polr0yvycZJkuPtdO8cg6oFEQXdJ56KyYaiqqi2lhhapNSJ/VUlPMyILdRAxkkHJ7NcZQLVAhgNCp07hNGmlZW9BwrJyi47v3j8hJ2Vnx4lKTGEK8fJVKmNESmIkqCczHJyIsNRRhyE8G29O3i4uMwxmuJlaOdM2raM8M8XlmkT3QE1rDWek7+Dpxhzp9ENojelqgot461Ehpt80XnsQXbc4HMZ6t0jkpw8F3v3F1JQWLZdbLRMb9tYkUPmEuSQ2UdF/oUul86uX/LIzStbFxKiccX52YSrGZVew7earOI+uhlTh0sahCCejhCCrBwX6RmOcqPrlJKe6WDfgUIcTl8vWlNhyIA44mPNxMfJLYCNDSmIJUhBBLvDy+GjRUHn0rxewXfbTuJylf26JLaxMra/CyX3l0qvoUUmYWxR/Vzb9YHLrXM8zV6p76IuBGknHBw+WkxcjJmO7X2+lKrqi78YG2NGVRXyC9z8tDuPPj0jCZXRuxskUhBLaO6CmJPr4tgJO6KCHtGBw0UcOFTkf68ocEG/GPp31XEf2wpU/FVSNJPPvUZrXGG28vJdpJ1w4C5nr3YpQvhW5LUz5k6NBpUW8RZ+3JXDjz/noakKXTqF0a93FDFRcnjdkJCCWEJzFURdF6SdsJMdxAm5TFsB237IprDQg8WsMvaSViTGg+vQtwg9SC9K6AhnAYolAgBjy15oEYm1fQvnBK9XkJHpICun7N7oquD2CL7dnIl+hqa2TQihe5dwOiTbMBjk8Lq+kf32Zkyx3cOx4/ZyYwWeiapA987hHEsrZszIVoRawHV4YzliKPCk78KbewhDfE8MrfuihjfelVdNU2jV0kpkpInjJ+zBw5xVwLHjxWXEEODw0WIOHy3GZFLp1N5Gt84RtG5paXDzqs0F2UMsoTn1EHVdlLiYlO3tZGY5cTi8JLYpG8RAUaBFvJW4GDNCCNzHtqEXBVlEEQJP5h68Wb/7DzS2HYo1ZeI52Y1S1wghyM1zcyLDUW7IszMptns5fLSYtPTyIwOVEhFupFvncLp1Dpd7ps8xUhBLaC6CWF6v0OXS+X1/AemZThQF+veJxmY7NYAwGlQS24T4FwM8mXvwZO8ve4EzxRBQLJEoplBC+t+GFtaybm6sHvB4dNIzHeTkuqo8jHa5dI6mFXP0uL3cnS+nk9DaSrfO4XTuGIZBk0PqukYKYglNXRDL6xV6vYLDR4s4dDRwSGcL0Ti/Xwyq4vO9S2wT4nch8eYdxX3ix7IXCSaGpjAUUwiWHhMwtupTR3dXvzgcvvzP1Ymk4/X6knCdSHdUOn9rNqncMrm9nGM8B8g5xGZAfoFveHd60ichIC3dt0c3mCtNYbGXjEwHKd0jaBF3ak7Lk70PT2YQ95ogYojB6hPDblc2WTEEsJRE0iko9D3nyiLpgG9OsmW8hZbxFhxOLycynJxIt1NsL3ts545hUgzPEVIQmzDFdg8n0h0UnbEAkJXtYu+BgjJJ4EvRNIWuncLo3zeaiHCfi4zQvb5FkvyyeZWDiqFmRDGH+cSwzXm1dk8NmTCbbwtgbp6b9IzK3XRKsZg1khNDSE4MKdlX7SA90+EfUsdEmcgvcBNmM5RZbNm0LQuTUSGpbSjRkTLlwdkih8wlNKUhs9Pp5USmo8xuk8IiD3v3F1Y4RGvV0kL3zhF06mDzx/oTHhfu49vQ7TllD/C6cZ/4ET3/2KkyRUOxRmPtMb7ZiOGZ6LrgZLaTrGxXlRdeyh7vIjvHRddOYYDvhyoizEhEhInQEA2HU2fJG/v8Ux1hNgPJbUNJTgwloU0IJrlLptpIQSyhKQiiy62TGSSPiNPlZf/BItLSHeUeGx1lolN7G+2SbMTHmv09Dd2Zj/vYNoS7bJQX4cjDfXQrwl10WqmCYo3C2nNisxXD0/F6Bdm5Lk5mOWskjOVhNKikZzr4YVdu0HpNVWjdykpyYojsPVYDKYglNGZBLCh0k53jyxZ35qd58EgxBw8XBvWBA9/iSccOYbRuaSWhdQjW0xLBewsz8KTtQOhnDK2FwJtzEE/GLl/W+dNQLBFYU/4kxfAMdP2UMLrdtSOMW7/PDhpcIhhhNgNtE0Jo0yqENq2sQYffEimIfhqbIHq9gpw835DKWYFj9b4DhRw6WjYyjcmk0D7JRuuWVuLjLGf0CgvwnPwFvTB4oAZP5q94T5ZdWFGs0Vj7/hljfPca3lXTR9d9n9vJLGeF8RerQmaWb0iele3EWc1zRUYY+fO1yVIUz0AuqjQidF1QWOQhv8CXlN3rFRQV+97bbEbCbWU/zqS2oRw/YcddMkGvqdA2MZS2JT6Fp/cKhceB5+RvePOOUtG+ZC2yLd6c/eAtnYtU0GI6EHL+NFRrVG3fdpNCVRViosxER5ooKPT4e/Y1IS7GTFyMGQijsNhDdraLk9lO8vLLjhTOJDSk/B6i0+nFZFKbpVhKQWzgeLw6BYUe8vNdZJ50kZvvpqDATV6Bm4JCt38o3DYhhHBb2Yx1Bk2hXbKN3/YW0LqllfZJoVitmj8Ki6IoCK8bb/Z+vDkHEKJylxHFaMXYqi/uo5tBUTG1H4mlx/hzngKgMXN6Dmmny0t2joucXFeN87PYQgzYQnzDYo9XkJPrqrD3aDGrHDpahNWiYbVomM2afxFmw7cZ7D9URGy0megoEzFRJqKjzMREmwixak1aKOWQuYSGMmT2egXFdl/P4Xi6ncxMJ/mFbvLz3f5eXjAiI4z06xW8d6YLsNs92EKNREeZiI81o2mgF51EL0jzJX3Sy85FCVchiqn8tKDenIOY2g7G2Lpv9W9UUgZdF+Tlu8nOdVV7r3RFFBZ7yMlxkZvnJifPhccj6JMSSXRkYNQhVVUwm1TWf5OO3RF8CG4xq0RHlQhltInoyKYllFIQS6gPQdR1gdOl43B4Kbb7XGKOp9vJz/dUOeBCKZoKwwbHlxeomsgIE/GxJoyeHLwFaegFJ4IHZcC3euzJ/AW98ASm5KEo1ugybdSQWIyt+6Focq9tXeB0ecnLd5OX5672d6Eyioo9WC1a0KyDDqeXjVuyqn1Os0khJtpCn56RdOoQVhtm1gtyyFzHCCHweAQut47LpeNwenE6dZwub5lJ9cwsJxmZzmpfw2JWCQ8z4vEKjIZTX3LFa8dmtBNnc2DUCxCH83CVI4IIHeHIxZO1F70gzV/szvwFU9vBAU21yGQM8d1QFOnnVleYTRrxsRrxsRYcDi+5+b4eXm2sUIeGlP/fPr+gZj1Tp0tw/ISd2GgTJrOKQVMwGFQ0TSn5t8KW7VmEhhqJjDASHWUmItzY4HwlpSDWkFKh83gEHq+OxyvwegQej47bI3C7dVxuLwUFHvIKfAsf+QVuEILz+pTtcQGEhxk4lha0yo/B4Jt7iggzEGYzEGETmFUXiteB4spFtTtQvcWEGIoIs+iYFA2KoMx/I68L3Z6Lbs9GL85GOLLhTPcaQBRlIoqzUEJiUBQNQ3w3tMikmj00SY2wWDRaWqy0iLPgcOgUFLopKPSUm/flbIiPNTNkQAwFhR6Kir0UFbl9f4s95bpunY7RqAZNP+H1Cn7ak1+m3GRSCA3x7fAJsxkICzUSFm4gJsqEQVNRVQVN8825qqqCqvqG9qpy2r8ryK9dXRqEIBYWFjJr1ixWrFhBdnY2Xbt25W9/+xt/+tOfKj02IyOD6dOns2rVKoqLi+nduzfz589n5MiRtWZfqZ+fx+sTPK9XBJ38drl18gvcFOS7yCtwU1joxu3RURCUrtqq6OguC5rqK1OE7qsTOtFmFzayUNBR0NFUHVuIQngI2KyCUKuOxSBAuFF0N4rXBXmnvqWaqhASYiA0zFAStblkkUMIhNuObs9G2LPRi7MQzrJfzqBoJnSPA1N4AobYzihG61k9S0nNURQFq1XDatWIj/NF2yko9FBQ6KawyFPjBZkzMZs0zNEasWf8btsdXoqKPBQWeykuDi6U5fU+7eXs73a5BC6Xb0HpdEZcGIdajTlJVVWwWjTaJ5c/510VGoQgjh8/nq1bt7JgwQI6d+7MW2+9xXXXXYeu61x//fXlHud0Ohk5ciS5ubksWrSI+Ph4XnzxRcaMGcPatWsZPnx4rdjndhTh3v8VCgIDAgO++T+Hw4PD4cXu9OCw637xMwAxJS+CjAi0tBCs1rIrshYB3aIcWMwaFouG2aTC6d8JAZT++ArdF6FB6JgMghALWEPMqKayc3q6PQv3oW+rd9OqES2mI6aEARha9kQ1N3zfzOaGwaASFWkiKtKEEMInWMVeios9FNu9tbozBvCvSMfGnCoTwhftp6jYg93hxWAILmLBglaUh8WsVksMwff/sTZWQ+p9UeWzzz7jsssu84tgKZdeeik///wzhw8fRtOCu3O89NJL3HnnnWzcuJFBgwYB4PF46N27Nzabjc2bN1fZjooWVfJz8sj88Qv/+6JiL0ePne7s7HuECr7enlryV0GU9PbEaXWCsOgIomKCTzx7cw/7hq5CL9kF4vvry3Wil6QF1TEqHsyaE4vqxKj6hk5aXDcMsZ2DnNSD87fPqMi3sBTFGIoakYixVR+MrXqhhsRWeoyk4SGEb8GuuNg39LU7vBU68Nc1eQVu0k44KLZ7SmwpX6yjIk30TYms9jVCQw20T2rkPcSVK1dis9mYOHFiQPlNN93E9ddfz+bNmxk8eHC5x3bp0sUvhgAGg4EbbriBmTNncuzYMdq0aXPWNp4ZmNNsVIhlv0/olFLxKx9FKZnr0BQ0VcFotgDBBVE4CxCi7NyQAhgVFxaDC4vmRFOCfLnL8yHUDCiWCIQj94yTqiiWSFRrNEpINKo1Gi00Di26A1p46wruSNLQURTFN9Iwa0SXeGN5vQKH0yeODscpkTwXXaKIMCMRp0X/1nVfj9Zu91Jc8tdu91Js9xISZPRUFarbqwxGvQvirl276NatGwZDoCm9evXy15cniLt27WLo0KFlykuP/fnnn8sVRKfTidN5akU3P7/8ObUzY9EZjCoG1VsmiTsAJeKnqaCpKmrJhPDpiqkoFXwDFQWEr7lB8WBU3ZhUNybViVbRcUBFs96qNRqvuxjVGo1aIn6KJRJUDcVkQwtrhRrWCtXceF0mJBWjaQqhIYaAeT5d93lAOJ06TqcXp+vU38pyUp8NqlrWlto459lS74KYlZVF+/bty5RHR0f76ys6trRddY998sknmTt3bpVsPDOtJICiqQiP17/ypWkqmgqKGih+QRFlhUtVFIwmFavZicGbj1Fxo1YmgFU4bymG+G4YWvSk1FFRMYaihZeKoJwfbK6o6qmeJJzqwQW4i7l13C7d/2+XS8fjEcE7BPWIWgsePPUuiECFHu6Veb/X9NgZM2bwwAMP+N/n5+eTmBg8RWbpMr9+2hfAajGCrgRdNDntyJLe4am/qqpgMhswhRgwGlQMhlP+WgBepQ3o8aBo+LqXvr+KooGqlvw1gKIGvkehPK9sRVFRQmJRLBGolgjfMFn2BCUVoCgKRqOC0agSGqReCJ+nhdstcHt0v7uZx+1zQavMI6MuaBI9xJiYmKA9uezsbICgPcDaONZsNmM2VzEDXInP0+lueoboUl88FRTFJzqaiqZqaAYNzaCiaT4HVa30pSq+HmQFaFHJVbMpmJmqhmIM8YXuN1pRzeG+BE9mm3SiltQqiqKU/JiDlYrn/HTdJ4qeEnE8831pWWm5rgu8ukDXwauLkgXFymkSPcSUlBTefvttPB5PwDziTz/9BEDPnj0rPLa03elU5djqYrVq6PopR1BNjS75e8phtNx9c9VG8QVKUDVf71I1+kLy+/8aQDP5ts0ZrChGC4rBKrfRSRokpf9njDX8egrhE0ddF+hCIHTQRamI4v9rsZy9Ita7283q1asZN24c77zzDtdee62/fOzYsezcubNCt5uXX36ZO+64g02bNjFgwADA53bTp08fbDYbmzZtqrIdFbndCN2LHpBLpHR4GrBSUrbutDa+IbNCaY8SRSv5e2o4fWp4LHtzEkm9IBoAo0aNElFRUWLx4sVi/fr1Ytq0aQIQb775pr/N1KlThaZp4uDBg/4yh8MhevToIRITE8W///1vsWbNGnHVVVcJg8EgNmzYUC0b8vLyBCDy8vJq7b4kEknjot6HzAAffPABjzzyCKmpqf6te2+//XbA1j2v14vXG+jqYjabWbduHdOnT+fuu++muLiYPn36sHr16lrbpSKRSJoP9T5kbig0lHiIEomk/pCTVRKJRFKCFESJRCIpoUHMITYESmcOKtrCJ5FIGjdhYWEVbtiQglhCQUEBQLm7VSQSSeOnsjUCuahSgq7rHD9+vNJfEDi1ze/IkSPNegFGPodTyGfho6E/B9lDrCKqqpKQkFCtY8LDwxvkh36ukc/hFPJZ+Gisz0EuqkgkEkkJUhAlEomkBCmINcBsNjN79uyqR8tposjncAr5LHw09ucgF1UkEomkBNlDlEgkkhKkIEokEkkJUhAlEomkBCmI1aCwsJD77ruP1q1bY7FY6NOnD++88059m1WnFBQUMH36dC699FLi4uJQFIU5c+YEbbtjxw4uueQSbDYbkZGRjB8/nv37959bg+uA9evXM3XqVLp27UpoaCht2rThiiuuYPv27WXaNtVnUMoPP/zAZZddRtu2bbFarURHRzNo0CDefPPNMm0b47OQglgNxo8fz/Lly5k9ezarV6+mf//+XHfddbz11lv1bVqdkZWVxeLFi3E6nVx55ZXltvvll18YMWIELpeLFStW8Prrr/Pbb78xdOhQMjMzz53BdcDLL7/MwYMHuffee/nss89YtGgRGRkZDBw4kPXr1/vbNeVnUEpubi6JiYk88cQTfPbZZ/zrX/8iOTmZyZMnM3/+fH+7Rvss6i82bePi008/FYB46623AspHjRolWrduLTweTz1ZVrfoui50XRdCCJGZmSkAMXv27DLtJk6cKGJjYwMijh88eFAYjUYxffr0c2VunZCenl6mrKCgQLRo0UKMHDnSX9aUn0FlDBgwQCQmJvrfN9ZnIXuIVWTlypXYbDYmTpwYUH7TTTdx/PhxNm/eXE+W1S2KolS6t9vj8bBq1SomTJgQsF0rKSmJiy66iJUrV9a1mXVKfHx8mTKbzUb37t05cuQI0PSfQWXExsb6k8Q15mchBbGK7Nq1i27dugVkBgTo1auXv765sm/fPux2u/9ZnE6vXr3Yu3cvDoejHiyrO/Ly8tixYwc9evQAmt8z0HUdj8dDZmYmL730Ep9//jkPP/ww0LifhRTEKpKVlRU0z3NpWbD80M2F0nsv7/kIIcjJyTnXZtUpd955J0VFRTzyyCNA83sGd9xxB0ajkfj4eO6//36ee+45brvtNqBxPwsZ7aYaVDR0rGxY2RxoLs/n0Ucf5d///jfPP/885513XkBdc3kGM2fO5JZbbiEjI4NPPvmEu+66i6KiIv7617/62zTGZyEFsYrExMQE7QVmZ2cDwX8NmwsxMTFA8F5ydnY2iqIQGRl5jq2qG+bOncv8+fN5/PHHueuuu/zlzekZALRt25a2bdsCMG7cOABmzJjBjTfe2KifhRwyV5GUlBT27NmDx+MJKP/pp58A6NmzZ32Y1SDo0KEDVqvV/yxO56effqJjx45YLJZ6sKx2mTt3LnPmzGHOnDnMnDkzoK65PIPyuOCCC/B4POzfv79xP4t6XuVuNHz22WcCEO+8805A+ZgxY5q0283pVOR2c80114j4+HiRn5/vLzt06JAwmUzi4YcfPodW1g2PPfaYAMSsWbPKbdPUn0FFTJ48WaiqKjIyMoQQjfdZSEGsBqNGjRJRUVFi8eLFYv369WLatGkCEG+++WZ9m1anfPbZZ+K9994Tr7/+ugDExIkTxXvvvSfee+89UVRUJIQQYs+ePcJms4lhw4aJzz77THzwwQeiZ8+eonXr1v7/JI2VZ555RgBizJgx4rvvvivzKqUpP4NSpk2bJh588EHx7rvvig0bNoj3339fXHvttQIQDz30kL9dY30WUhCrQUFBgbjnnntEy5YthclkEr169RJvv/12fZtV5yQlJQkg6OvAgQP+dtu2bRMjR44UISEhIjw8XFx55ZVi79699Wd4LTF8+PBy7//MQVZTfQalvP7662Lo0KEiNjZWGAwGERkZKYYPHy7eeOONMm0b47OQ8RAlEomkBLmoIpFIJCVIQZRIJJISpCBKJBJJCVIQJRKJpAQpiBKJRFKCFESJRCIpQQqiRCKRlCAFUdIoKA1UW9lrw4YNTJkyheTk5Po2WdIIkY7ZkkbBpk2bAt7PmzePL7/8MiCnCUD37t3JzMwkPz+fvn37nksTJU0AGf5L0igYOHBgwPu4uDhUVS1TDgSErZdIqoMcMkuaHMGGzIqicNddd7F06VK6dOmC1Wrl/PPPZ9OmTQghePrpp2nXrh02m42LL76YvXv3ljnv2rVrGTlyJOHh4YSEhDBkyBDWrVt3ju5Kci6QgihpNqxatYolS5awYMEC3n77bQoKCrjssst48MEH+fbbb3nhhRdYvHgxu3fvZsKECZw+m/Tmm29y6aWXEh4ezvLly1mxYgXR0dGMHj1aimITQg6ZJc0Gp9PJF198QWhoKODrNV555ZV8+eWX7Nixwx/WPjMzk/vuu49du3aRkpJCcXEx9957L5dffnlAxrhx48bRr18/Zs6c2WSzLjY3ZA9R0my46KKL/GII0K1bNwDGjh0bkOOjtPzQoUMAbNy4kezsbG688UY8Ho//pes6Y8aMYevWrRQVFZ3DO5HUFbKHKGk2nJn3xmQyVVhemiozPT0dgKuvvrrcc2dnZweIraRxIgVRIqmE2NhYAJ5//vmgq9oALVq0OJcmSeoIKYgSSSUMGTKEyMhIdu/eHZBpT9L0kIIokVSCzWbj+eef58YbbyQ7O5urr76a+Ph4MjMz+fHHH8nMzOTll1+ubzMltYAURImkCtxwww20bduWp556ittuu42CggLi4+Pp06cPU6ZMqW/zJLWE3LonkUgkJUi3G4lEIilBCqJEIpGUIAVRIpFISpCCKJFIJCVIQZRIJJISpCBKJBJJCVIQJRKJpAQpiBKJRFKCFESJRCIpQQqiRCKRlCAFUSKRSEqQgiiRSCQl/D+FnFS3CCorfQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(3.2, 3))\n", "ax.plot(time, cell_number[:, 0], '--', lw=3, label='Stem\\nsimulated', c='#9098d9')\n", "ax.plot(time, cell_number[:, 1], '--', lw=3, label='Non-stem\\nsimulated', c='#ed9e44', alpha=0.7)\n", "ax.plot(range(T), n_stem.sum(0), label='Stem\\nestimated', c='#9098d9', lw=6, alpha=0.4)\n", "ax.plot(range(T), n_nonstem.sum(0), label='Non-stem\\nestimated', c='#ed9e44', lw=6, alpha=0.4)\n", "ax.legend(loc=2,fontsize=10)\n", "ax.ticklabel_format (style='sci', scilimits= (-1,2), axis='y')\n", "ax.spines['right'].set_visible(False)\n", "ax.spines['top'].set_visible(False)\n", "ax.set_xlabel('Time')\n", "ax.set_ylabel('Cell number')" ] }, { "cell_type": "markdown", "id": "a3b2dcac-4848-4268-807e-b46df2f05821", "metadata": {}, "source": [ "Accurate fitting of LP and LR distance distributions can illustrate the accuracy of the inference" ] }, { "cell_type": "code", "execution_count": 9, "id": "f83ece39-da76-43e4-b49d-b8821a085e0e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "lp_Power_divergenceResult(statistic=0.13853071748300083, pvalue=0.999638361551232)\n", "lr_Power_divergenceResult(statistic=0.008769196264595613, pvalue=0.9997821707208558)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAEoCAYAAAAHaLa7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB/MklEQVR4nO3dd1iTV/sH8G/IBAJhK0MZooIKolatWgdF3K0VtRVt3zra2j3ft7W1rtbZ+murdrr1tWq1YlsXCq66UFwoioqLIaLMhAAJJDm/P3jz1MiGQEJyf64rV8uTZ9wJ4XjnOec+h8cYYyCEEEIIIWbJxtQBEEIIIYSQ6lGyRgghhBBixihZI4QQQggxY5SsEUIIIYSYMUrWCCGEEELMGCVrhBBCCCFmjJI1QgghhBAzRskaIYQQQogZo2TtfxhjUCgUoDmCCSHmhtonQqwbJWv/U1RUBJlMhqKiIlOHQgghBqh9IsS6UbJGCCGEEGLGKFkjhBBCCDFjlKwRQgghhJgxStYIIYQQQswYJWuEEEIIIWaMkjVCCCGEEDMmMHUAhOhptVqUl5ebOgzSQgiFQvD5fFOHQZpATEwMFixYgGvXriEoKAgzZ85EVFSUyeJhjEGr1UKj0ZgsBtKyGLt9omStnsytEbEEjDFkZ2ejsLDQ1KGQFsbJyQmtW7cGj8czdSjESGJiYjB27FgMGTIEX375Jfbv34+xY8dix44dzd7WMsZQWFiInJwcaLXaZr02afmM2T7xGE2JDQBQKBSQyWSQy+VwdHSscp+YmBi88MILePbZZ9G7d28cPHgQBw4cMEkjYknu37+PwsJCeHh4wM7Ojv7hJbVijKGkpAQPHz6Ek5MTPD09TR1Sk6pL+2QpevToAZlMBplMhhdffBFRUVEYOnQo8vPzcfbs2WaNRd82OTo6wtHREQKBgNonUqumaJ8oWfufujSGPXr0gKOjI/r06YN79+7hp59+wnPPPWeSRsRSaLVa3LhxAx4eHnB1dTV1OKSFycvLw8OHD9GhQweL7hK1pmTN3t4e/fv3x/79+9GpUydcuXIF33zzDWbNmoXi4uJmi0Or1SI1NRVubm5wc3NrtusSy2HM9okKDOrh2rVrCA8Px+bNm5GQkIDbt29j2LBhSElJMXVoLZZ+jJqdnZ2JIyEtkf5zQ2MdLUdQUBCOHz8OAEhJSYFCoUBsbCyCg4ObNY7y8nIwxmBvb9+s1yWWw5jtEyVr9RAUFIQTJ06gTZs2AIAjR46YpBGxRNS1QBqCPjeWZ9KkSdwdNMYYIiIiEBcXh5kzZ5okHvqMkYYy5meHCgzqYebMmRg7diw6deoEAFi5ciUuX76MmJgYE0dGCCGWITc31+Dn7OxsxMTEYMyYMSaKiBDToztr9RAVFYUdO3ZAP8xPrVbj999/p0aEEEKMgDGGrVu3AgC6d+8OAOjbty+1scTqUbJWT1FRUThy5AhsbCreug4dOpg4ImKu1q9fDx6Pxz0kEglat26N8PBwLFq0CA8fPmzQea9evYq5c+fi7t27xg2YEBNLTEzEnTt3YG9vj7lz53LbiPFR+9SyULLWAA4ODvDx8QEAxMXFmTgaYu7WrVuHU6dOIS4uDj/88APCwsKwZMkSBAcHIz4+vt7nu3r1KubNm0eNIbE4W7ZsAQA8++yz6N+/PwDgzp07lbpGifFQ+9QyULLWQPo7agkJCSaOhJi7Ll264Mknn0T//v0xduxYfPvtt7h06RLs7e0RFRWFBw8emDpEQkxOq9Xit99+AwBMmDABTk5OXDtLUyM1HWqfWgZK1hooNDQUAHDp0iXodDoTR0MeFxMTgx49esDe3h49evQwuyKQtm3b4v/+7/9QVFSEX375BUDFP0gTJkyAn58fbG1t4efnh+joaKSlpXHHrV+/HuPHjwcAhIeHc10Y69evB1Bxp3f06NHw8fGBRCJBYGAgpk+fTncmiNk7fvw47t+/DycnJwwdOhQA8MQTTwCwvK5Qap+ofaovqgZtoKCgIIjFYhQXF+PChQvo0aOHqUOyGPrZnxvqzz//xKRJkxAREYFZs2YhPj4eY8eOxa+//orRo0c36JxNsbLCiBEjwOfz8ffffwMA7t69i44dO2LChAlwcXHB/fv38dNPP6Fnz564evUq3NzcMHLkSCxcuBCfffYZfvjhB24Qdrt27QAAt27dQp8+ffDKK69AJpPh7t27+Oabb/DUU0/h8uXLEAqFRn0NhBiLvrAgKioKYrEYANCzZ09s3rzZrJI1c2ufmmrVF2qfzAwjjDHG5HI5A8Dkcnmt+5aUlLCEhAT23HPPsQ4dOrAFCxY0Q4SWqbS0lF29epWVlpZy25RKJQNgVg+lUlnv17Zu3ToGgCUmJla7T6tWrVhwcHCVz2k0GqZUKpm9vT1btmwZt3379u0MADt8+HCN19fpdKy8vJylpaUxAOzPP/+s92swd1V9fixRfdqnlqisrIy5uroyAOzAgQPc9uPHjzMAzNPTs9ljqu6zZW7tU0PaJsaofWoOxmyfqBu0gXg8HoKCggCA++ZBSH2xR1Z7UyqV+OSTTxAYGAiBQACBQACpVIri4uI6r5Lx8OFDvP7662jTpg0EAgGEQiF8fX0BgFbaIGbr4MGDyMvLg7u7O8LDwwFUjGELCwsDn8/H/fv3ce/ePRNHaX2ofTIf1A3aCCEhIYiJiUFycjLUajV36540jp2dHZRKZYOPf+qpp+Dq6oo//vgDPB4PjDGMHj0aBQUFOHbsWINjMrbi4mLk5eUhJCQEADBx4kQcPHgQs2bNQs+ePeHo6Agej4cRI0agtLS01vPpdDoMGTIEWVlZmDVrFkJCQmBvbw+dTocnn3yyTucglSmVSnz++efYtm0b8vPzERQUhBkzZmDChAn1Os/nn3+OBQsWoHPnzkhOTm6iaFsmfRfo+PHjIRAIUFJSgvT0dNjb26Nz5864dOkSEhMT4e3tbeJIza99aqql+qh9Mi+UrDVCQEAAHBwcUFxcjFOnTmHQoEGmDski8Hi8Rq3HN2vWLIwdOxZRUVEYNmwYYmNjcejQIcTExJjVOn979uyBVqvFoEGDIJfLsXv3bsyZMwczZszg9lGr1cjPz6/T+ZKTk5GUlIT169fj5Zdf5rbfvHnT6LFbk6ioKCQmJmLx4sXo0KEDNm/ejOjoaOh0OkycOLFO57h48SKWLl2KVq1aNXG0LY9KpcLOnTsBANHR0QCAnJwcFBQUICsrC0888QSXrD333HMmjLQCtU/UPpkCdYM2gkgkQocOHSAWi2m+NTOiX2kiPz8fs2bNQn5+vtktV5Oeno5///vfkMlkmD59OvcN+/G7s6tXr4ZWqzXYpt/n8W+i+kHGj59DX81F6m/v3r2Ii4vDjz/+iOnTpyM8PByrVq1CZGQk/vOf/1T63VRFo9FgypQpmD59Ojd0gvxj3759UCgU8PHxQd++fVFUVISioiIUFBSgtLQUYWFhACynIpTaJ2qfGsIs7qw1tpvhzz//xDfffIMLFy5Aq9XCz88P7733Hl577bUmjVsoFCIoKAjnzp3DoUOHmvRapH6ioqIQFRVl6jAAVHyj1Gg00Gg0ePjwIY4dO4Z169aBz+dj586dcHd3BwAMGDAAX3/9Ndzc3ODn54ejR49izZo1cHJyMjhfly5dAFSsTevg4ACJRAJ/f38EBQWhXbt2mDFjBhhjcHFxwa5du+iLRCPs3LkTUqmUm45Ab8qUKZg4cSJOnz6Nvn371niOxYsXIz8/HwsWLMCoUaOaMtwWSd8F+sILL4DH4yEvLw9yuRxqtRo8Hg/t27cHUDF1BGPMIhZWp/aJ2qf6Mos7a1FRUdiwYQPmzJmDffv2oWfPnoiOjsbmzZtrPXbx4sWIiopCly5dsG3bNvz111948803UVZW1uRx83g8rj//ypUrkMvlTX5N0vJMmTIFffr0QUREBN544w1cuHABn3zyCa5du8YNpgaAzZs3Izw8HB9//DGioqJw9uxZxMXFQSaTGZzP398f3333HZKSkjBo0CD07NkTu3btglAoxK5du9ChQwdMnz4d0dHRePjwYYNmIScVkpOTERwcDIHA8Hutfp7F2saeXb16FfPnz8dPP/0EqVTaZHG2VEqlErt27QJQMRGu/q6aQqGAs7MzAMDDwwNisRgFBQW4deuWKcO1SNQ+tRCNridtpD179jAAbPPmzQbbIyMjmZeXF9NoNNUee/bsWWZjY8OWLFnS6DjqO3XH6dOn2eXLl9mpU6dYv379mKurK/vjjz8aHYe1sZapF0jTaOrPT/v27dnQoUMrbc/KymIA2MKFC6s9VqvVst69e7Po6Ghu28CBA1nnzp1rva5KpWJyuZx7ZGRkWOTUHZs3b2YAWGBgINNqtezWrVvs2LFjbP/+/ezChQvsr7/+YkeOHGG9e/eu8t+JpkRtE2ksi5q6o6ZuhqysLJw+fbraY7///nuIxWK88847TR1mtYRCIQIDAyGRSOgbAiEWqKZut5qe++abb5Camorvvvuu3tdctGgRZDIZ92jTpk29z9ES6NcCnTBhAhQKBZRKJRQKBVxdXSEWiyEUClFcXMxNOm4p49YIqS+TJ2uN6Wb4+++/ERwcjB07dqBjx47g8/nw8fHBjBkzmqUbFKgoMujYsSPEYjEla4RYGFdXV+Tl5VXarq+Ac3FxqfK49PR0zJ49G3PmzIFIJEJhYSEKCwuh0Wig0+lQWFhY41QFn376KeRyOffIyMgwzgsyIwUFBYiNjQVQMV5NP1ZNIBDA3t4eNjY2kEgkKC8vR+fOnQFQskasl8mTtby8vCobPP22qhpKvXv37iE1NRXvvvsu3n33XcTHx2Py5MlYunQppkyZUuN11Wo1FAqFwaMhBAIBgoKCIBQKkZqaShM3EmJBQkJCkJKSAo1GY7D98uXLAP4ZTP2427dvo7S0FO+99x6cnZ25x4kTJ5CSkgJnZ2d8+umn1V5XLBbD0dHR4GFpdu7cifLycoSEhMDb2xtKpRJFRUXw8PDg7lja2tqCMQY/Pz8AwPnz5yv9LgixBiZP1oCGdzPodDoUFRXhxx9/xFtvvYXw8HDMnz8f77zzDjZv3lzj/C3G6mawsbGBi4sLtzDtwYMHG3QeQoj5GTNmDJRKJXbs2GGwfcOGDfDy8kLv3r2rPC4sLAyHDx+u9OjatSv8/Pxw+PBhvP32283xEszWo1WgeXl5KCgogEgkgq2tLbePSCSCjY0NHB0d4eDggJKSEprpnlglkydrDe1m0B8LAEOHDjXYPnz4cAAV38KqY8xuBrFYTOPWCLFAw4cPR2RkJN544w2sWrUKhw8fxmuvvYbY2Fh89dVX4PP5AIBp06ZBIBAgLS0NAODk5IRBgwZVejg5OcHe3h6DBg1CYGCgKV+aST148ID7YjtixAgolUqUlJRw00To6cetFRUV0bg1YtVMnqw1tJsB+Gdc2+PY/9Yzs7Gp/uUZs5tBKBSiffv23J019sh6aoSQli0mJgYvvfQSZs+ejWHDhuH06dPYsmULJk2axO2j1Wqh1Wrpb7+Ofv/9d+h0OvTq1QtSqRSFhYWQSCSQSCQG+/H5fIjFYpSXl1vc5LiE1IfJk7WGdjMAwNixYwFUzID9qL1798LGxgY9e/Y0fsBVEIlEaNeuHaRSKbKysnDt2rVmuS4hpOlJpVIsW7YM9+/fh1qtRlJSUqUJu9evX28wtqo6R44coXVBYbgWaFFREUpKSuDm5lblvra2ttDpdOjYsSMAStaIdTL5CgaPdjMoFAoEBgZiy5YtiI2NxaZNmwy6GTZs2IBbt27B19cXQMX0Hr/88gvefPNN5ObmolOnToiPj8cPP/yAN998k9uvqQkEAkgkErRv3x43b95EfHw8goODm+XahBDSkmRkZOD48eMQCATo27cvCgsLYWdnV2kZIj2xWAwejwcfHx8AwKVLl6BWq6vdnxBLZPI7a0DDuxmEQiHi4uIwYcIELFy4ECNGjMDOnTuxePFiLFu2rNnit7Gxga2tLdq1a0dFBoQQUoPffvsNADB48GBIJBKoVKpq76oBFT0XAoEAAoEAbm5uKC8vR1JSUnOFS4hZMItkrTHdDC4uLvj555+RnZ2NsrIyXL9+Hf/+979rHK/WFB4tMjh8+DCVlxNCSBW2bt0KPp+PyMhIFBQUwN7eHkKhsNr9BQIBRCIRysvLuaEt1BVKrI1ZJGuWQCgUwtvbG25ublAoFDh79qypQyImtn79evB4vGofR44cqdN5srKyMHfuXFy8eLHSc3PnzjXZwtZXr17F3LlzcffuXaOf25SvizSd1NRUnDt3Dh4eHggLC4Nara7xrhpQMX2TnZ0dtFotV3BGyVrjUfvUcKZ4XSYfs2YpRCIRRCIRQkJCcP36dRw8eBBPPvmkqcMiZmDdunUICgqqtL1Tp051Oj4rKwvz5s2Dn58fVxGn98orr2DYsGHGCLPerl69innz5mHQoEG1DqwnBKjoAhUKhejfvz90Oh2kUmml1Wuqoh+fFhAQAICSNWOi9qlloGTNSAQCAfh8Pjp06AA+n4/4+HjMnDnT1GERM9ClSxc88cQTTXJuHx8fbuA1IeaMMYYtW7bA1dUVvXr1QllZGTw9Pet0rH7cWqtWrQAAKSkpKCoqgoODQ1OGbBWofWoZqBvUSB4vMjh58iSKi4tNHRZpAbZv347evXtDJpPBzs4OAQEBmDp1KoCKqR7043SmTJnCdVHMnTsXQNW34/38/DBq1Cjs3r0b3bp1g62tLYKDg7F7924AFd0fwcHBsLe3R69evSp12Z89exYTJkyAn58fbG1t4efnh+joaG7CV/05xo8fDwAIDw/n4lq/fj23T3x8PCIiIuDo6Ag7Ozv069evyuKbPXv2ICwsDGKxGP7+/li6dGnj3lBilpKTk3Hz5k24uLigXbt2kMlkdbqrBlQka0KhECKRCL6+vmCM1TjpOTEeap/Mo32iZM2IxGIxnJ2d4efnh7KyMhw/ftzUIbVYOp3OrB6NodVqodFoDB5arRYAcOrUKbzwwgsICAjA1q1bsWfPHsyePZsrUOnevTvWrVsHAPj8889x6tQpnDp1Cq+88kqN10xKSsKnn36KTz75BDExMZDJZIiKisKcOXOwevVqLFy4EL/++ivkcjlGjRplsKj43bt30bFjR3z33XfYv38/lixZgvv376Nnz57Izc0FAIwcORILFy4EAPzwww9cXCNHjgQAbNq0CUOGDIGjoyM2bNiAbdu2wcXFBUOHDjVoEA8ePIjRo0fDwcEBW7duxddff41t27Zxr5lYjq1bt8LNzQ1du3aFSCSCs7NznY/l8XiwtbVFeXk5unfvDsC0XaGmbo+M1TYB1D61lPaJukGNSCgUwsbGBj169MCVK1dw8ODBSkthkdrpdDqkpqaaOgwD7du3b3CFcVVjF/l8PjQaDU6ePAnGGH7++WfIZDLu+cmTJwMAHB0duUHV7dq1q/M4yLy8PCQkJMDb2xsA4OXlhbCwMKxatQo3b96EnZ0dgIp/BJ977jnEx8fjmWeeAQCMGzcO48aN486l1WoxatQotGrVCps3b8a7774Ld3d3tG/fHkDF2JZH4yopKcF7772HUaNGYefOndz2ESNGoHv37vjss89w+vRpAMDMmTPRqlUrxMXFcbPXDx061CLGmJB/MMawY8cOSKVShISEwNnZmZtDs670n4/g4GDs3LnTZMmaubVPjWmbAGqf9My9faI7a0akH1ehH6xJ64QSANi4cSMSExMNHvrGQN+F8Pzzz2Pbtm24d++eUa4ZFhbGNYQAuEmaBw0axDWEj25/tAtBqVTik08+QWBgIDe/lVQqRXFxcZ0W0T558iTy8/Px8ssvG3xb1+l0GDZsGBITE1FcXIzi4mIkJiYiKirKYJkhBwcHrmEmliExMRFKpRIikQjBwcFwcnKq9zlEIhH4fD73uaYiA+Og9qlltE90Z82I9B8cPz8/CAQCXLhwAbm5ubWWphNDNjY23Lcic9GYb67BwcHVDuAdMGAA/vjjDyxfvhz/+te/oFar0blzZ8ycORPR0dENvqaLi4vBzyKRqMbtKpWK2zZx4kQcPHgQs2bNQs+ePeHo6Agej4cRI0YYdEdU58GDBwBg8O33cfn5+eDxeNDpdGjdunWl56vaRlqu3377Dfb29ujSpQs8PT0b9PckFou5IgM+n487d+6YpH01t/apsXOKUvtUmTm2T5SsGZGNjQ0kEgns7OwQEhKCCxcu4PDhw9xAR1J3zT2psSmNHj0ao0ePhlqtRkJCAhYtWoSJEyfCz88Pffr0adZY5HI5du/ejTlz5mDGjBncdrVajfz8/DqdQ/+P54oVK6rtFmnVqhXKy8vB4/GQnZ1d6fmqtpGWSafT4fDhwwCAkJAQODo6Nug8+iKusrIyBAUF4cqVKzh79qxJpoag9onap+ZmPZ+4ZiIWi8EY4z7E1BVK6kosFmPgwIFYsmQJAODChQvcdgB1+tbYWDweD4yxSusurl69mht0/Gi8VcXVr18/ODk54erVq3jiiSeqfIhEIq7aKyYmxuCbc1FREXbt2tVEr5A0t0OHDqG4uBhisRgDBgxoVKIjkUjAGENoaCgA0OTjzYjapwqmap/ozpqRCYVC8Hg8btAlJWskOTm5yuXH2rVrhxUrViAzMxMRERHw8fFBYWEhli1bBqFQiIEDB3L72dra4tdff0VwcDCkUim8vLzg5eVl9FgdHR0xYMAAfP3113Bzc4Ofnx+OHj2KNWvWVBpnpP+Mr1y5Eg4ODpBIJPD394erqytWrFiBl19+Gfn5+Rg3bhw8PDyQk5ODpKQk5OTk4KeffgIAfPnllxg2bBgiIyPx0UcfQavVYsmSJbC3t6/zN2Vivhhj+OuvvwBUDIRvbJelWCyGjY0NTY5rRNQ+tYz2ie6sGZm+yKBNmzYQCAS4ffs27ty5Y+qwiAlNmTIFffr0qfT4888/0bt3b2RnZ+OTTz7BkCFD8Nprr8HW1haHDh1C586dAQB2dnZYu3Yt8vLyMGTIEPTs2RMrV65ssng3b96M8PBwfPzxx4iKisLZs2cRFxdnUA0GAP7+/vjuu++QlJSEQYMGoWfPntw3zhdffBGHDx+GUqnE9OnTMXjwYLz33ns4f/48IiIiuHNERkbijz/+gEKhwAsvvIAPP/wQY8eO5eZxIi1bQUEBTp06BZ1Ox8131Rj69tXLyws8Ho+SNSOg9qlltE88xhhr9quaIYVCAZlMBrlcXuuYitLSUly+fBl2dnaVFiDW6XRIT0+HRqPB4sWLcfz4caxatarWeWeslUqlwp07d+Dv729QcUNIXVjL56c+7ZO5YIxh+/btmDVrFsrKyrBr164aF2yvq3v37qGgoACTJ0+GUqlEZmamQWWhsVjLZ4s0HWN+hujOmpHpiww0Gg2eeuopANQVSgixPgqFAgcPHoRWq0XPnj2NkqgBgK2tLfh8Prd2Jd1dI9aAkrUmoC8y0M+0fejQIaPMNE0IIS2BTqfD/fv3cfr0aeTn5+PZZ5812rn10zno57OkZI1YA0rWmoC+yMDHxwdSqRQ5OTm4fPmyqcMihJBmIZfLcezYMRQVFUEsFqNHjx5GO7dYLIZQKISvry8AStaIdaBkrQnoB8EqlUquYoa6Qgkh1kCn0yEvLw9Hjx5FXl4ehg8fbtR5yQQCAUQiEby8vCAWi3H27FnQ0Gti6ShZawL6lQzKysowYMAAAJSsEUKsQ0FBAXJycpCYmAi5XI4RI0YY9fw8Hg92dnbw8PCATCZDQUEBbt26ZdRrEGJuKFlrAo8WGejXVvv7779RVlZm4sjMF30zJg1BnxvzotVqkZ+fj+PHjyM7Oxtt27bl5rsyJv18a81RZECfMdJQxvzsULLWRPRFBp6envDw8EBJSQkSEhJMHZbZEQgq5mWualJGQmqj/9zoP0fEtAoKCqBUKnHs2DEoFAqMHDmy0XOrVUU/1CQwMBBA0yRrfD4fAFBeXm70cxPrYMz2iZK1JqIvMigqKuIm2aOu0Mr4fD74fD4UCoWpQyEtkEKh4D5DxLQ0Gg3y8vKQlZWFY8eOAYDRu0D1hEIhhEIhvL29IRAImiRZEwqFEIvFkMvldHeNNIgx2yf6OtpE9N/8CgsLERERgS1btuDgwYP44osvTB2aWeHxePDw8MD9+/chFothb2/fJN/EiWVhjKG4uBgKhQKenp70mTED+fn5KC4uxqlTp1BYWIj27dujffv2TXIt/aLuXl5esLW1xfnz56HRaIx+h9XNzQ337t1DZmYmZDIZ9yWckJo0RftEyVoTEQqFXJFB3759AQCnT5+GQqFoMTOQNxeZTIbS0lLk5uYiJyfH1OGQFoLH48HJyanSMjOk+ZWXlyM/Px+FhYU4cuQIAGDkyJFNek2JRAI3Nze4urri7t27SElJQUhIiFGvoW+rc3Nzce/ePaOem1g2Y7dPlKw1ER6PB4lEgsLCQshkMgQGBuLmzZs4evQonnnmGVOHZ1Z4PB43to/Gh5C6EgqF1P1pJvLz86FUKvHw4UOcOnUKQNN1gerpey+CgoJw9+5dJCYmGj1ZAyoSNkdHR5SXl0Or1Rr9/MQyGbt9omStCemLDIqKijB48GDcvHkTBw8epGStGjT2iJCWp6ysjLurdu7cOeh0OoSEhKBt27ZNel2RSMRNjmtjY4PExMQmXWBbP06OEFOgAoMmpB/foFAoMHjwYABUZEAIsSx5eXkoLi5GeXk5Dh06BKDp76oBFV/ubG1t4e3tDVtbW1rJgFg0Staa0KNFBoMGDQKPx8OVK1dw//59U4dGCCGNplarUVhYiIKCApSWluL8+fPg8XgYPnx4s1xfIpFwydqlS5egVqub5bqENDdK1prQo0UGYrHYYGF3Qghp6XJzc1FUVASNRsPNI9mjRw+0bt26Wa4vEong4uLCjXdNSkpqlusS0twoWWtC+iIDjUZDXaGEEIuiUqkgl8tRUFAAmUyGffv2AWj6KtBHicViiEQiBAYGgsfjUVcosViUrDWxx4sMgIpkjSZZJIS0ZLm5uVAqldDpdJDL5UhOTgafz8eQIUOaLQb9ou4+Pj4Qi8WUrBGLRclaE3u0yKBfv34Qi8XIzMzEjRs3TB0aIYQ0SElJCeRyOfLz8+Hs7Iz9+/cDAJ588km4uro2ayx2dnbw9vaGnZ0dJWvEYlGy1sQeLTKQSCTo168fAODgwYMmjowQQhomNzeXWyJOJpNh7969AJq3C1RPLBajTZs2sLW1RUpKCoqKipo9BkKaGiVrTezRIoPi4mIat0ZIC6NUKvH+++/Dy8sLEokEYWFh2Lp1a63HxcfHIzIyEl5eXhCLxfDw8MDTTz/NJTYtVWlpKYqKiiCXy+Hi4oKbN28iNTUVQqGQa9+ak77IwNPTE4wxnD9/vtljIKSpUbLWxB4vMtAv6n748GGaDZuQFiAqKgobNmzAnDlzsG/fPvTs2RPR0dHYvHlzjcfl5eWhc+fO+Pbbb3HgwAH88ssvEAqFGDlyJDZt2tRM0RufSqVCaWkpdDodHB0dsWfPHgDAgAEDTLKUnn6yWl9fX4hEIuoKJRaJVjBoBo8WGfTo0QMymQyFhYU4f/48evbsaerwCCHV2Lt3L+Li4rB582ZER0cDAMLDw5GWlob//Oc/eOGFF6pddeOFF17ACy+8YLBt1KhR8Pf3x8qVK/Hiiy82efxNQa1Wo6ysjJvNX3+nsDkmwq0Kj8fjxq3R5LjEUtGdtWagLzKQy+Xg8/l4+umnAVBXKCHmbufOnZBKpRg/frzB9ilTpiArKwunT5+u1/mEQiGcnJwgELTc78kqlQoqlQpisRjJycnIyMiAra0twsPDTRaTWCyGj48PJWvEYlGy1gz0RQZyuRyMMa4rlJI1QsxbcnIygoODKyVXoaGh3PO10el00Gg0yMrKwpw5c3Djxg189NFHTRJvU9PpdFCpVFCr1bC1teW6QMPDw2FnZ2eyuMRiMfz8/GBra4s7d+4gLy/PZLEQ0hQoWWsG1RUZnDhxAqWlpSaOjhBSnby8PLi4uFTart9Wl6RgxIgREAqF8Pb2xnfffYfffvut1qpJtVoNhUJh8DAHZWVlKCsrg06ng1gsNslEuFURCoVwcHCAl5cXBAIBzp49a9J4CDE2StaaweNFBh06dICPjw/UajVOnDhh6vAIITXg8XgNek5vxYoVOHPmDP78808MHToUL7zwArZs2VLjMYsWLYJMJuMebdq0qXfcTUGlUqGsrAwAkJSUhAcPHsDBwQH9+/c3aVw2Njaws7OjrlBisShZayaPFhnweDzqCiWkBXB1da3y7ll+fj4AVHnX7XHt27dHz5498eyzz2Lbtm2IiIjAW2+9BZ1OV+0xn376KeRyOffIyMho+IswIrVajfLycgiFQu6uWmRkJEQikYkjq1jUnZI1YqkoWWsmIpGIKzIAwHWF0uS4hJivkJAQpKSkQKPRGGy/fPkyAKBLly71PmevXr1QUFCAnJycavcRi8VwdHQ0eJgDtVoNtVoNGxsbxMbGAjBdFejjRCIR2rZtS8kasUiUrDUT/bi1x4sMzp07x31LJ4SYlzFjxkCpVGLHjh0G2zds2AAvLy/07t27XudjjOHo0aNwcnJq9mWZGosxhtLSUpw+fRqvvfYaCgoKwOfzzWbFgEeTtQcPHuDevXumDokQo2m59eMtzONFBp6enujcuTOuXLmCw4cPY+zYsaYOkRDymOHDhyMyMhJvvPEGFAoFAgMDsWXLFsTGxmLTpk3cHGvTpk3Dhg0bcOvWLfj6+gIARo8eja5duyIsLAyurq7IysrC+vXrcfToUfzwww8tbvqOsrIy/P3331i9ejWUSiUAwN3dHR988EGzL+BeFT6fD5lMhlatWiE9PR2JiYnw9vY2aUyEGAvdWWsmjxYZ6LtCadwaIeYvJiYGL730EmbPno1hw4bh9OnT2LJlCyZNmsTto9VqodVqwRjjtvXr1w+xsbF45ZVXEBERgXfeeQc8Hg+7d+/Gm2++aYqX0ihqtRo7duyAv78/SkpKAABff/01+vbti19++cXE0VWwtbWlcWvEIlGy1oweLTIAaNwaIS2BVCrFsmXLcP/+fajVaiQlJWHChAkG+6xfvx6MMfj5+XHbPv74Y5w5cwb5+fnQaDTIzc1FbGysyae5aCiVSoV79+7B2dkZxcXF8PDwQI8ePdC/f3/cvn3b1OEBqOgK1S/qTskasSSUrDUjfZGBfs6kgQMHgs/nIzU1FWlpaSaOjhBCqqdWq9G6dWtcvHgRANC7d2/weDwcO3YMAQEBpg3uf0QiETc57tmzZw3udBLSklGy1oweLzJwdHREr169ANDdNUKIeSstLcXgwYO5L5YqlQrTpk3DyZMn8frrr5s4ugpCoRB+fn4QCoUoLS3FrVu3TB0SIUZByVozerzIAKCuUEKI+SsrK4NarUbnzp25iYD//vtvyOVyrFixApGRkSaO8B9SqRSenp7UFUosSoOTNf0s1qTuqioy0Cdr8fHxdMueECOh9sm49JPhZmdnQ6VSwdbWFufOncOOHTvMKlEDKsYG07g1YmkanKx5e3vj008/RXp6ujHjsXiPFxk8+eSTsLOzw8OHD+u0KDQhpHbUPhmXSqWCRqPhuhW7dOlitlOPiMVi+Pr6UrJGLEqDk7VnnnkGy5cvR7t27TBmzJhGdeMplUq8//778PLygkQiQVhYGLZu3Vrv83z++efg8XgNmlW8uTxeZCASiTBgwAAA1BVKiLEYs30i/6xccOPGDQBAWFiYaQOqgVAohL+/P/h8Pi5fvlxp9QlCWqIGJ2tr165FZmYmFixYgKSkJAwZMgTBwcH4/vvv6z2jdVRUFDZs2IA5c+Zg37596NmzJ6Kjo7F58+Y6n+PixYtYunQpWrVqVd+X0qweLzIADLtCCSGNZ8z2iVTcWSstLcWVK1cAAN26dTNxRNXj8Xjw8/ODSCQCYwwpKSmmDomQRmtUgYGzszM+/vhj3Lp1Czt37kSbNm3w3nvvwdvbG2+//TauXbtW6zn27t2LuLg4/Pjjj5g+fTrCw8OxatUqREZG4j//+Q+0Wm2t59BoNJgyZQqmT5+OoKCgxrykJldVkYF+ctwjR46gvLzclOERYjGM0T6RivZVrVZDLpdzd9a6du1q4qhqZmdnB29vb9jZ2VFXKLEIRqkG5fF4ePbZZ7FkyRIMHDgQSqUSP/74Izp37oyxY8fi4cOH1R67c+dOSKVSjB8/3mD7lClTkJWVhdOnT9d6/cWLFyM/Px8LFixo9GtpalUVGYSGhsLNzQ3FxcV1er2EkLprTPtEKu6qlZWV4e7du2CMoW3btma/rikt6k4sTaOTNY1Ggy1btuCpp57CE088gdu3b2PJkiW4e/cuvvvuOxw7dgz/+te/qj0+OTkZwcHBlQarhoaGcs/X5OrVq5g/fz5++uknSKXSxr6cZvF4kYGNjQ13d43G1hBiPI1tn8g/laA3b94EYN7j1fREIhH8/f0hFApx9uxZU4dDSKM1uJzn3r17+OWXX7Bq1So8ePAA/fv3x7Zt2zBmzBjY2FTkgO+88w68vb3x4osvVnuevLy8Kme/dnFx4Z6vjk6nw9SpUxEVFYURI0bUK379gFk9/YD/5vB4kQFQMW7tt99+Q3x8PObMmdNssRBiiYzVPpF/7qxdv34dQMtI1mxsbBAYGAgAuHXrFtRqNcRisYmjIqThGpys+fn5QSAQYMKECXjvvfeq/QMOCAioddC/fpLF+j73zTffIDU1FX/99VedYn7UokWLMG/evHofZwyPFxnweDzuzlpCQgKUSmWLuUtIiDkyZvtk7dRqNYqLi7nigpaQrAFAmzZtYGdnB4FAgKSkJG61GEJaogZ3g86ZMwfp6elYt25djX+8YWFhuHPnTrXPu7q6Vnn3LD8/H8A/d9gel56ejtmzZ2POnDkQiUQoLCxEYWEhNBoNdDodCgsLUVpaWu11P/30U8jlcu6RkZFR7b7GVlWRgb+/PwICAqDRaPD33383WyyEWCJjtU/WTqvVorS0FJmZmcjPz4ednR06dOhg6rDqhCbHJZakwcla27Ztue6Ex+Xn52Pjxo11Ok9ISAhSUlIqzYVz+fJlAKh2zrTbt2+jtLQU7733HpydnbnHiRMnkJKSAmdnZ3z66afVXlcsFsPR0dHg0VyqKjIAaAoPQozFWO2TtVOr1dxkuDqdzqwnw32cfnJcsViMM2fOmDocQhqlwcnalClTql0k986dO5gyZUqdzjNmzBgolUrs2LHDYPuGDRvg5eWF3r17V3lcWFgYDh8+XOnRtWtX+Pn54fDhw3j77bfr96Ka0eNFBsA/U3hQskZI4xirfbJ2+vFqLam4QI/P5yMgIAA8Hg9JSUmmDoeQRmnwV6Sa1rFUqVTg8/l1Os/w4cMRGRmJN954AwqFAoGBgdiyZQtiY2OxadMm7jzTpk3Dhg0bcOvWLfj6+sLJyQmDBg2qdD4nJydoNJoqnzMnVRUZPP300wAq7io+ePCAxtIQ0kDGap+snVqtRllZGTexrDlPhluVjh07AqgoOCkqKoKDg4OJIyKkYeqVrKWnp+Pu3bvczxcuXIBKpTLYp7S0FCtXrkTbtm3rfN6YmBjMnDkTs2fPRn5+PoKCgrBlyxZMmDCB20er1UKr1VrMYudVFRm4ubmhW7duuHDhAg4dOoTo6GhTh0lIi9FU7ZM1U6lUyMvL48b1taQ7awDQunVryGQylJSU4Pz58xg4cKCpQyKkQeqVrK1btw7z5s0Dj8cDj8fDm2++WWkffTK1bNmyOp9XKpVi2bJlNR6zfv16rF+/vtZzHTlypM7XNaXHiwz01Z8RERG4cOEC4uPjKVkjpB6aqn2yVjqdDiqVCjdu3IBKpYKvry+cnZ1NHVa96CfHLSgowOnTpylZIy1WvZK1559/Hl26dAFjDM8//zwWLlyI9u3bG+wjFovRpUsX+Pn5GTNOi6MvMigsLIRcLueStcGDB2Pp0qWIj4/n7rgRQmpH7ZNx6btAb9++Da1W2+LuqgGAQCCAr68vLl++jHPnzpk6HEIarF7JWnBwMIKDgwFUfIsdNWqU2S87Ys6qKjJ46qmnIBKJkJ6ejlu3bnETOxJCakbtk3GpVCqUl5dz64G2tPFqQMWXYn0bevXqVRNHQ0jDNbga9OWXX6aGsJGqKjKwt7dHnz59AFBVKCENRe1T46nVaqhUKm7B+5Z4Zw0Al8Dn5eXVuCIOIeasXnfWvvjiC7zyyivw8vLCF198UeO+PB4Ps2bNalRwlq6qIgOgoiv06NGjiI+Px+uvv27iKAlpGah9Mi6VSoXbt29DLpfDzs6uUpdyS+Hi4gI3NzdoNBokJiZi2LBhpg6JkHrjsXqUV9rY2CAhIQG9evWqdsJJ7sQ8HrRabaMDbC4KhQIymQxyubzWCXJLS0tx+fJl2NnZQSgUNviajDGkpaVBo9Fg4MCB3Li1hIQE9OnTBy4uLnj48CFNM0BIHVD7ZDw6nQ43btzApk2bsGTJEjzxxBN1KvAyR4wxfPHFFzh37hyio6Mxd+5cU4dESL3V686aTqer8v9Jw1RXZPDEE0/A0dER+fn5uHjxInr06GHiSAkxf9Q+GU9ZWRlXXKDRaFpsFyhQ0c76+fnh3LlzuHjxoqnDIaRBGjxmjRhHVUUGAoGAm9SXxq0RQpqbWq02KC5oyckaAG4909TUVBNHQkjDNDhZU6lUBgPjAWDbtm2YMWMGJRj1UFWRAfDPOqEHDx40RViEtGjUPjWOSqVCfn4+0tPTAQBdu3Y1cUSNExQUBBsbG5SUlODevXumDoeQemtwsvbSSy/h3Xff5X5evnw5JkyYgK+++gpDhw7F3r17jRKgpXu8yEBPn6wdO3as0izshJCaUfvUOGq1GikpKVCpVPDz82txk+E+TiaToXXr1hCJRDhx4oSpwyGk3hqcrJ05c8agqmb58uV48cUXUVhYiKioKCxdutQoAVq6x1cy0AsKCoKnpydUKhVOnjxpwggJaXmofWo4xhhKS0uRmpoKtVrd4rtAgYrikzZt2gAAEhMTTRwNIfXX4GQtJycH3t7eAIA7d+7g9u3beOedd+Do6Ihp06YhOTnZaEFaMn2RgUajgVwuN9hOXaGENAy1Tw2nLy64desWysrKLCJZA4CAgAAAwJUrV0wcCSH11+Bkzc7Ojksujh07BqlUiieeeAIAIJFIoFQqjROhFaiqyAD4pyuUxtgQUj/UPjWcWq2GWq3mBuO3xJULqhIUFAQAuHv3LuoxYxUhZqHByVpISAh++OEHXL58GT/++CPCw8O5SV3T09PRunVrowVp6aorMoiIiAAAnD17FoWFhSaIjJCWidqnhtNPhqtQKGBvb28xS94FBQWBz+dDo9FwqzIQ0lI0OFmbNWsWjh49irCwMCQlJeHjjz/mntuzZw+6d+9ulACtQXVFBt7e3vD29oZOp0Pr1q3Ro0cPxMTEmDBSQloGap8aTq1W49q1a1CpVOjatavFTMpta2sLLy8v8Hg8KjIgLU69JsV91NNPP42UlBScO3cOYWFh3HgA/XOWMs6hOTxeZKCfHDcmJoYrM3/iiSdgb2+PsWPHYseOHYiKijJlyISYNWqfGk6lUuHGjRsWU1zwqLZt2yIjIwMXLlwwdSiE1EuDkzUA8PX1ha+vb6Xt06dPb8xprU51KxksWLAAYWFhuHjxInJzc3Hs2DEMHToUCxcupGSNkFpQ+1R/5eXlUKvVuHPnjkUma4GBgThx4gQ32S8hLUWjkjUAePjwIdLS0lBaWlrpuQEDBjT29FajqiKDa9eu4bPPPsPly5dx/fp1ZGRkYNiwYbQANSF1RO1T/ahUKuTm5nJ39Fv6ZLiP69KlCwAgKysLarUaYrHYxBERUjcNTtbu37+Pl156CYcPH670HGOsxS2UbGpVFRkEBQXh77//Ru/evXHy5Ens27cPsbGxCA4ONmGkhJg/Y7ZPSqUSn3/+ObZt24b8/HwEBQVhxowZmDBhQo3HxcTEYPv27UhMTMS9e/fQqlUr9OvXD3PnzkX79u0b9Lqa2qPj1QICAiCTyUwdklG1a9cOQqEQ5eXlOHfuHPr27WvqkAipkwYna2+//TYuXLiAJUuWIDQ0lL6hNNLjRQY8Hg8zZ87E2LFj0a5dOwDA7Nmz8fDhQyoyIKQWxmyfoqKikJiYiMWLF6NDhw7YvHkzoqOjodPpMHHixGqPW7JkCVq3bo2ZM2ciICAAGRkZWLhwIbp3746EhAR07ty5wTE1FZVKhevXr0OtVqNfv36mDsfoBAIBPD09kZ6ejjNnzlCyRlqMBidrR48exdKlSzFlyhRjxmO1qioyiIqKwo4dO/D5558DAHJzc7Ft2zaMGTPGxNESYt6M1T7t3bsXcXFxXIIGAOHh4UhLS8N//vMfvPDCC9VWS+7atQseHh4G255++mn4+fnh22+/xerVqxsVW1NQqVRITU2FSqWyuPFqegEBAUhPT8fly5dNHQohddbgqTt4PB63fAdpvOpWMoiKikJycjLc3d2h0+kqNf6EkMqM1T7t3LkTUqkU48ePN9g+ZcoUZGVl4fTp09UeW9XfqpeXF3x8fJCRkdHo2IxNo9GgpKQEd+/etcjiAr2OHTsCAG7fvm3iSAipuwYna+PHj8fu3buNGYvVq24lAxsbGwwdOhQAsG/fPlOERkiLYqz2KTk5GcHBwRAIDDshQkNDuefr4/bt20hLSzPLLlC1Wo1bt26huLjYoibDfVxoaCgYY3jw4AHy8vJMHQ4hddLgbtDnn38er776KnQ6HZ555hm4urpW2ocmnqyf6lYyAIDhw4dj06ZNiI2NxeLFi00QHSEth7Hap7y8PIM52vRcXFy45+tKo9Fg2rRpkEql+OCDD2rcV7/kk15VbYKxqVQqpKSkQKVSITQ0FDY2Df4ub9Z8fHwgFAqh0Whw+vRpjBgxwtQhEVKrRk2KCwDff/89fvjhB4PnqBq0YaoqMtAbMmQIeDwekpKSkJWVBS8vLxNGSoh5M2b79OjfYX2ee/ya06ZNw7Fjx7Bjx45au2gXLVqEefPm1encxqJWq7niAkvtAgX+6SK/c+cOzp8/T8kaaREanKytW7fOmHEQVL+SAQC4ubmhZ8+eOHPmDPbv30+FHYTUwFjtk6ura5V3z/Lz8wH8c4etJowxvPLKK9i0aRM2bNiA0aNH13rMp59+ig8//JD7WaFQNPkYYZVKhZs3b0KtVlvM4u3VCQwMxJ07d3D9+nVTh0JInTQ4WXv55ZeNGQfBP0UGBQUFBisZ6A0bNgxnzpxBbGwsJWuE1MBY7VNISAi2bNkCjUZjMG5NX0mon2S1OvpEbd26dVizZg1efPHFOl1XLBY363RIWq0W2dnZyM3N5bpBLVmnTp0QFxeH9PR0qFQqSCQSU4dESI2MMijh+vXrOHHiBIqLi41xOqumb6AfLzIAKsatAcCBAweg0WiaNS5CWqrGtE9jxoyBUqnEjh07DLZv2LABXl5e6N27d7XHMsbw6quvYt26dfjll1/M+guWWq3G1atXUV5eDj8/P4ubDPdxoaGhKCsrw4MHD8yyMpeQxzUqWdu4cSN8fHzQqVMnDBgwgLul/Pzzz2PVqlVGCdDa1FRk0LNnTzg7O6OwsBBnzpwxQXSEtBzGaJ+GDx+OyMhIvPHGG1i1ahUOHz6M1157DbGxsfjqq6+4OdamTZsGgUCAtLQ07th3330Xa9aswZQpUxASEoKEhATuYW4LiatUKly7ds0qukABwN3dHQKBAIwxDBo0CD169KDJxolZa3Cytn37dkyePBndu3fH999/D8YY91z37t2xbds2owRobR4vMngUn8/HkCFDANAUHoTUxJjtU0xMDF566SXMnj0bw4YNw+nTp7FlyxZMmjSJ20er1UKr1RpcZ9euXQCAtWvXok+fPgYPc5vYWq1W48aNGxY9Ge6jDhw4gAcPHgAAevfuDTc3N4wdO5YSNmK2GpysLVq0CFOmTMFff/2F1157zeC54OBgXL16tdHBWaNHiwyUSmWl5/VdobGxsc0dGiEthjHbJ6lUimXLluH+/ftQq9VISkqqtC7o+vXrwRiDn58ft+3u3btgjFX5uHv3bmNentEVFxfjzp07Fl8Jqvfzzz/D1dWVqwz+7bffEBkZiYULF5o6NEKq1OBkLSUlpdqFjF1cXGiywQbSFxmUl5dX2RWqnxz37NmzePjwYXOHR0iLQO1T3el0OqSkpKCsrAwCgYBbi9iS3blzB4MGDYJSqcS1a9dw584dDBs2DCkpKaYOjZAqNThZs7OzM1gW6VH37t2Ds7Nzg4OydjUVGbRu3ZobU3LgwIFmjYuQloLap7rTFxdoNBp06dLFYifDfZS/vz9SU1Ph6OgIjUaD/fv3IzY2FsHBwaYOjZAqNfivsl+/fpXGguitX78egwYNakxcVq2mIgOgYgoPgLpCCakOtU91p1arraq4AABef/11rkJYq9Vi27ZtOH78OGbOnGnq0AipUoOTtdmzZyMhIQG9evXC8uXLwePxEBMTg2eeeQZ///03fegboaYiA+CfcWv79++HTqdr7vAIMXvUPtWdSqVCamqq1RQXABUrwixfvhwSiQRKpRLFxcVYtGiR2RV+EKLX4GTtiSeewL59+6BUKvHRRx+BMYaFCxfixo0b2Lt3b62TRZLq1VZk8OSTT8LR0RG5ubk4d+6cCSIkxLxR+1R39+/fR05ODsrKytC1a1dTh9NshgwZgt27d8Pd3R1AxSoR9OWXmKtGDU4IDw9HSkoKUlNTcfz4cVy7dg3Xr1+nLoZGqq3IQCgUYvDgwQBoCg9CqkPtU+10Oh2SkpIAAJ6ennBwcDBxRM3v6aefhkajQWJiInJyckwdDiFValCylpOTg/nz52PIkCHo0qULRo8ejblz52LHjh1UZWUkNRUZADSFByHVofap7srKynD16lVotVqEhISYOhyTGD58OIqKipCamoobN26YOhxCqlTvtUEPHjyIsWPHQqFQgM/nw83NDYwxXL9+HfHx8Vi6dCl27tyJAQMGNEW8VkNfZFBdRZu+yOD06dPIz8+v04LShFg6ap/qR61W4/r161ZVXPA4f39/eHp6QqVSIT4+Hn369DFYB5YQc1CvO2s5OTl44YUXIJPJsG3bNsjlcty/fx/Z2dmQy+XYunUr7O3tMW7cOPoG20j6cWsKhaLKIgMfHx906dIFOp0OcXFxJoiQEPNC7VP9KZVK3Lp1y6qTNQAYPHgwysrKaP5KYrbqlaytWbMGWq0WJ06cwLhx42BnZ8c9Z2dnh+effx7Hjx9HeXk51qxZY/RgrUltRQbAP3fXaNwaIdQ+NcSVK1dQXl4OgUAAf39/U4djMsOHD4dCocCtW7doYlxiluqVrB04cABTp06Fj49Ptfu0bdsWU6ZMobFUjcTj8WBra1ttkQFgOG6NqpiItaP2qX4YY7h06RIAIDAw0Comw61OmzZt0LZtWzDGcOjQIZSVlZk6JEIM1OuvMyUlBU899VSt+/Xv35++nRiBSCQCUH2RQb9+/WBvb48HDx5wjS4h1orap/opLy9HSkoKGGNWW1zwqCFDhkClUuHcuXPIzs42dTiEGKhXslZYWAgPD49a9/Pw8EBhYWFDYyL/U1uRgVgsRkREBADqCiWE2qf6UalUuH79OlQqFbp3727qcExO3xV69+5dXL161dThEGKgXsmaWq2GUCisdT/9WCvSOLUVGQC09BQhetQ+1c+9e/e4yXCtZeWCmnh6eiIwMBCMMRw5cgSlpaWmDokQTr3rk69fv15rWfO1a9caHBD5x+NFBlVNWKlP1k6cOAG5XA6ZTNbcYRJiNqh9qjv96ifu7u6QSqUmjsY8DB06FBs2bOC6Qq256IKYl3ona5MnT651H8YYeDxeQ+Ihj9AXGeTn50OhUFSZrPn7+6Njx464fv06Dh48iKioKBNESoh5oPap7pKTkwEAHTt2NHEk5mPo0KH4/vvvkZ6ejsuXL1OyRsxGvZK1devWNVUcpBq1FRkAFXfXrl+/jn379lGyRqwWtU9192hxQWhoqKnDMRseHh4ICgpCfn4+/v77b4SHh1vlElzE/NQrWXv55ZebKg5SjdqKDICKgbHLli1DbGws3TUgVovap7orKirCrVu3UFZWhieeeMLU4ZiV4cOHY+XKlTh//jyys7MpWSNmwXon1mkh9OPW5HJ5tXOpDRgwABKJBJmZmbhy5UozR0gIaWkuXrzITYYbEBBg6nDMypAhQ1BcXIx79+5xi9wTYmqUrJk5oVAIoVCIsrIyPHjwoMp9bG1tER4eDoCqQgkhtbtw4QIAwNfXl+7EP8bFxQWdO3eGTqfD8ePHUVBQYOqQCDGPZE2pVOL999+Hl5cXJBIJwsLCsHXr1lqPi4mJQXR0NAIDA2Fraws/Pz9MmjQJqampzRB18+DxeHBxcUF5eTlSU1NRXl5e5X40hQchpK70d+CpuKBqI0aMgFKpxIULF6r9kkxIczKLZC0qKgobNmzAnDlzsG/fPvTs2RPR0dHYvHlzjcctWbIEJSUlmDlzJmJjYzF//nxcuHAB3bt3t6juQKlUCplMhsLCQty5c6fKffTJ2rFjx6pdS5QQQjQaDW7cuAEA6Nq1q4mjMU+DBw9GaWkpsrOzcf78eVrOj5hcvafuMLa9e/ciLi4OmzdvRnR0NAAgPDwcaWlp+M9//oMXXngBfD6/ymN37dpVacbyp59+Gn5+fvj222+xevXqJo+/OfB4PLi7u0OpVCItLQ2enp6VBr22b98eAQEBuH37Ng4fPoxnnnnGRNESQsxZWloacnJyoNFoaOWCajg5OSEsLAyZmZk4deoUIiMj4e7ubuqwiBUz+Z21nTt3QiqVYvz48Qbbp0yZgqysLJw+fbraY6taWsbLyws+Pj7IyMgweqymJBKJ4ObmBqVSiRs3blRa0YDH43F312jpKUJIdfRtqpubGxwdHU0cjfkaMWIEioqKqCuUmAWTJ2vJyckIDg6uNOu4fu4f/cSNdXX79m2kpaWhc+fORovRXDg5OcHW1hYPHjzA/fv3Kz0/fPhwABXJWnXLUxFCrNvFixcBVBQXkOpFRERApVIhJycHZ8+ehUajMXVIxIqZPFnLy8uDi4tLpe36bXl5eXU+l0ajwbRp0yCVSvHBBx/UuK9arYZCoTB4mDs+nw93d3eUlZXh5s2blYoNwsPDIRKJcPfuXW5MCiGEPEq/SHlwcLCJIzFvUqkUvXr1Qnl5Oc6cOYPc3FxTh0SsmMmTNQA1lo7XtaycMYZp06bh2LFj2LhxI9q0aVPj/osWLYJMJuMete1vLuzt7eHo6IjCwkLcvn270nMDBgwAQFWhhJDKVCoV7t69CwC0eHsdDB8+HAqFAhcvXkR2drapwyFWzOTJmqura5V3z/Lz8wGgyrtuj2OM4ZVXXsGmTZuwfv16jB49utZjPv30U8jlcu7RUsa46YsNACA9Pb3SMlQ0bo0QUp3ExESUl5dDJBKhffv2pg7H7A0aNAjl5eXIy8vD2bNnq506iZCmZvJkLSQkBCkpKZXGA1y+fBkA0KVLlxqP1ydq69atw+rVq/Hiiy/W6bpisRiOjo4Gj5aipmIDfbJ29OhRlJaWmipEQogZOnPmDADAx8cHQqHQxNGYP3t7e/Tt2xdqtRpnzpyhQgNiMiZP1saMGQOlUokdO3YYbN+wYQO8vLzQu3fvao9ljOHVV1/FunXr8Msvv2DKlClNHa7Z0BcbZGdnGxQbdOrUCW3atIFKpcLRo0dNGCEhxNzol0/y8/MzbSAtiL4rNCkpibpCicmYPFkbPnw4IiMj8cYbb2DVqlU4fPgwXnvtNcTGxuKrr77i5libNm0aBAIB0tLSuGPfffddrFmzBlOmTEFISAgSEhK4h345FUulLzYoLy83KDagKTwIIdW5du0aAFhktXxTGThwILRaLQoLC5GYmAiVSmXqkIgVMnmyBlQsG/XSSy9h9uzZGDZsGE6fPo0tW7Zg0qRJ3D5arRZardagy2/Xrl0AgLVr16JPnz4GjzFjxjT762hu9vb2kMlkKCgoMCg20E/hQUUGhBC9zMxMFBYWgsfjcVMjkdpJJBIMHDgQpaWlOHv2LHWFEpMwi2RNKpVi2bJluH//PtRqNZKSkjBhwgSDfdavXw/GmMHt+7t374IxVuVDX/FkyXg8Htzc3MDj8ZCWlga5XA6gYn4ggUCAGzduVKoYJYRYpxMnToDH46FVq1Z1Ktwi/9B3hV66dAlZWVmmDodYIbNI1kjDiUQiuLu7o6SkBKmpqdDpdHB0dES/fv0A0N01QkiFxMREAECbNm0qTUJOavbUU08BAORyOc6ePYvi4mITR0SsDSVrFsDJyQkSiQQPHjzgBsDSuDVCyKP0xQX+/v51nr+SVBCJRHj66adRUlKCc+fOUaEBaXaUrFkAGxsbrtggNTUVZWVlXLJ26NAhqNVqE0dISMulVCrx/vvvw8vLCxKJBGFhYdi6dWutx2VmZuL999/HwIED4eTkBB6Ph/Xr1zd9wFVQq9XckIhOnTqZJIaWbsSIEVAoFLh8+TIyMzNNHQ6xMpSsWQh9sYFcLsft27fRtWtXtG7dGiUlJTh27JipwyOkxYqKisKGDRswZ84c7Nu3Dz179kR0dDQ2b95c43E3b97Er7/+CpFIhBEjRjRTtFU7d+4c+Hw+7O3tERAQYNJYWqonn3wSAoEACoUCFy5c4MYIE9IcKFmzEPpiA6BiZQOFQsHdXaNxa4Q0zN69exEXF4cff/wR06dPR3h4OFatWoXIyEj85z//gVarrfbYAQMGICcnB3Fxcfjwww+bMerKTp06BR6Ph7Zt28LW1taksbRUQqEQgwcPRnFxMXWFkmZHyZoF0RcbFBcX48aNGzRujZBG2rlzJ6RSKcaPH2+wfcqUKcjKysLp06erPdbGxnya17NnzwKoKC4QiUQmjqbl0neFJicn4969ewZTSRHSlMynNSFGoV/Z4OHDhwgJCYGNjQ2uXr2K9PR0U4dGSIuTnJyM4ODgStWT+nnKkpOTTRFWvemX7wsICKDigkbo2bMnJBIJioqKcO7cOW4Na0KaGiVrFubRYoOcnBz07dsXALB//34TR0ZIy5OXl1flnGT6bXl5eU1yXbVaDYVCYfBoqIyMDMjlctjY2CA4ONiIUVofgUCAoUOHQqlU4vz58zRBLmk2lKxZIHt7ezg5OaGwsBADBgwAQF2hhDRUTXeimuou1aJFiyCTybhHmzZtGnyukydPQiKRwMvLC05OTsYL0krpu0KvXr2KzMxM6HQ6U4dErAAlaxZIX2xgY2MDPz8/2NraIj4+nls/lBBSN66urlXePdN3fzXVSgCffvop5HI598jIyGjwuRISEmBjY0PFBUbSvXt3ODo6cnfXcnJyTB0SsQKUrFkooVAINzc3uLq6IjAwEEVFRTh58qSpwyKkRQkJCUFKSgo0Go3Bdv0YsC5dujTJdcViMRwdHQ0eDXX+/HkAQNu2bam4wAhsbGwwbNgwFBUV4cKFC9QVSpoFJWsWzMnJCfb29ujcuTOcnZ1pCg9C6mnMmDFQKpXYsWOHwfYNGzbAy8sLvXv3NlFkdaNSqXDjxg0AQLt27cyqQrUl068VmpKSgoyMjErJPCHGRgvEWTAbGxt4eHigQ4cOOHv2LPbv349FixaZOixCWozhw4cjMjISb7zxBhQKBQIDA7FlyxbExsZi06ZN4PP5AIBp06Zhw4YNuHXrFnx9fbnjf//9dwDgVg84e/YspFIpAGDcuHFNHv/58+fB5/Ph4OAAb2/vJr+etejatStcXV25u2vdunWDl5eXqcMiFoySNQtnZ2eHHj164LfffkN2djaysrKoUSGkHmJiYjBz5kzMnj0b+fn5CAoKwpYtWzBhwgRuH61WC61WW2nercfnZ/vhhx/www8/AECzzNF16tQpSCQS+Pn5wc7OrsmvZy14PB6GDx+OnTt34uLFi8jOzqZ2lTQpuiduBQICAuDj4wMHBwfs2rXL1OEQ0qJIpVIsW7YM9+/fh1qtRlJSkkGiBgDr168HYwx+fn4G2xlj1T6aQ0JCAvh8Ptq0aUPJmpENHz4cRUVFuHbtGjIyMlBWVmbqkIgFo2TNCgiFQm5+paNHj1KpOSFWgDGGCxcuAAB8fX2puMDIOnfuDE9PTxQVFXF31whpKpSsWYm+fftCo9EgJSWFVjMgxApkZGSgoKAANjY2CAgIoOICI9N3hSoUCiQlJVFVKGlS9NdrJUJDQ6FSqVBSUoKDBw9CrVabOiRCSBPSj1ejLtCmM2LECBQVFeH69evIyMhAaWmpqUMiFoqSNSshEAjQtWtXlJSU4Pz587h165apQyKENKGTJ09CLBZTcUET6tChA/z8/KBQKKgrlDQpStasSP/+/fHw4UNuQGxhYaGpQyKENJGEhAQIhUK0adMGEonE1OFYpMe7QilZI02FkjUr8tRTT6GsrAyXL1/Gw4cPcePGDSo2IMQClZaW4tq1awAqigsoWWs6I0aMgFKpxI0bN5CZmYmioiJTh0QsECVrVsTDwwPBwcHIzc3FnTt3kJOTg3v37pk6LEKIkZ07dw58Ph+Ojo5wc3PjJu8lxhcQEID27dtTVyhpUpSsWZmnnnoKjDEkJSWhvLwct27domIDQizMqVOnIBaL0a5dO5qyoxlQVyhpapSsWZkBAwYAAA4fPgyZTIbCwkIqNiDEwugrQakLtHmMGDECxcXFSE1NRVZWFgoKCkwdErEwlKxZmbCwMNjb26OgoAAPHz4En8/n5mMihLR8jDEkJCRAJBLRtB3NpG3btujUqRMUCgUuXLhAc64Ro6NkzcoIhUL07dsXQEVpv7u7O0pKSqjYgBALkZaWhoKCAvD5fHh7e0MsFps6JKswYsQIKBQKXLp0Cffv36f2lBgVJWtWqH///gCAv//+GzKZDHZ2dsjNzUVmZqaJIyOENJa+C7Rt27aQSCQQCASmDskqDBs2DCUlJUhNTcX9+/eRl5dn6pCIBaFkzQrpk7VLly5BoVDAw8ODig0IsQAxMTF47733IBaLkZubi/Pnz5s6JKvh7e2NsLAwKBQKnD9/nrpCiVFRsmaFPD090b59e+h0Opw8eRJ2dnZwdnaGQqHAzZs3TR0eIaQBYmJiMHbsWKjVakgkEjg4OGDRokU4cOCAqUOzGvqu0MuXLyM7OxsajcbUIRELQcmalXrqqacAAMeOHQMAuLq6csUG+fn5pgyNENIACxYswODBg1FcXAyRSIRp06aha9eu+OWXX0wdmtUYOnQo1Go11xWam5tr6pCIhaBkzUrpp/A4fvw4dDodhEIh3N3dUVpaitTUVBocS0gLc+3aNXTs2BECgQDOzs5wcnJC3759cfv2bVOHZjVatWqFHj16cFWhNOcaMRZK1qxUjx49YGtri5ycHFy/fh0A4OjoCHt7e+Tk5FCxASEtTFBQEG7cuIH4+Hi89tpr4PF4OHnyJAICAkwdmlUZMWIEioqKcOnSJTx8+BDl5eWmDolYAErWrJRIJELv3r0BVFSFAoCNjQ08PDyg0Whw8+ZNZGRkUENDSAsxc+ZMxMXFYenSpcjIyMDSpUtx8uRJvP7666YOzaoMGTKEa0Pv3btHhQbEKChZs2L6rlD9uDUAsLW1hYeHB+RyOVJSUpCUlIQHDx5Aq9WaKkxCSB1ERUVhx44dKC4uxo4dO1BYWIgVK1YgMjLS1KFZFTc3N/Tu3ZvWCiVGRRPwWDH9FB4XL15EUVERHBwcAFQUG9jZ2SE7OxuZmZlQKBRwc3NDmzZtIJPJYGNDOT4h5mjMmDHo3Lkzbt++DVdXV8hkMlOHZJWGDx+OxMREXLp0Cbm5uVCr1TQ5MWkU+lfXirVp0wZ+fn7QaDQ4deqUwXO2trbw9/dHmzZtUFJSgrS0NFy6dAnXr1+HUqk0UcSEkJqUlZWhrKwMWq2W1gQ1If3dzFu3biEzMxNpaWlUtEUahZI1K6e/u/ZoV+ijpFIpAgMD4ebmhsLCQty+fRtJSUm4ffs2VCpVc4ZKCKmFjY0NZDIZrVxgYs7Ozujbty8KCgqQlJSE9PR03LhxgyYdJw1GyZqVe3TcGmOsyn14PB5cXFzQvn17ODg4cBWkSUlJyMzMpCIEQsyEUCiEi4sLpFKpqUOxesOHD4dSqcSvv/6Kd999F2+99RZWr16NvLy8attaQqpDyZqV69mzJ8RiMbKzs2tdvYDP56NVq1YIDAyESCTC/fv3kZKSwpWoUxECIYRU0I/tLSgoQJ8+fSAUCrF8+XJs2LABt2/fpi+5pF4oWbNyEokEvXr1AvDPFB61EQqF8PHxgZ+fH7RaLTIzM5GcnIwrV66gsLCQvjUSQqzehg0b4OLiAgBISEjA5MmT4e3tjV27diE1NRXJyckoKioycZSkpaBkjXDj1o4fP16v4/RFCN7e3lSEQAghj7hz5w5GjhwJe3t7XL16FdHR0Xj48CHOnz8PxhgyMjK4oSTUK0FqQ8ka4dYJPXv2LIqLi+t9vIODA1eEUFBQgFu3blERAiHEqvn7++P27dvYt28fxo4dCwC4cuUKlEolDhw4wK0Wk5KSgmvXrqG0tNTEERNzRskagb+/P3x8fFBeXo7Tp0836ByPFiFIpVLk5OTgxo0buHTpEu7du0fjMwghVuX111/HiRMnMGPGDAQGBiI0NBQAoNPp8N1332Hy5MlIS0tDaWkp7t69i6SkJDx8+JCGkZAqUbJGwOPxGtwV+jg+n4/WrVujXbt2EAqFyMrKwtWrV3H58mUqQiCEWI0hQ4Zg+fLlkMvlWL58ObRaLZYvX46vv/4aHh4eyMzMxAcffIAlS5YgNzcX2dnZuHLlCm7evImysjJTh0/MDI9RGg8AUCgUkMlkkMvlcHR0rHHf0tJSXL58GXZ2dhAKhc0UYdM6dOgQ3nzzTfj4+CAuLg48Hs8o5y0tLcX9+/dRVlYGBwcHuLu7cyshGOsahFi6+rRPSqUSycnJcHR0BJ/Pb6YISX0UFxdj1apVWLt2LcrKysDn8zFp0iQMGzYMTk5O8PDwQGBgIK1AQTh0Z40AAHr37g2hUIjMzEzcvXvXaOe1tbVFQEAAvLy8UFxcjLt37+Ly5cu4ceNGg8bHEUJIS2dvb4/3338fe/fuRWRkJLRaLTZu3Ih33nkHf//9N9LS0pCUlIS0tDTqjSAAKFkj/2Nvb48nnngCQPWrGTSGo6Mj2rdvD1dXV+Tn5+PmzZu4ePEi7ty5Q0UIhBCr5OPjgxUrVmD9+vVo37498vLysHjxYsyfPx+nT5/G9evXceXKFfpiSyhZI//QV4U2RbIGVIyNc3V1Rfv27blKqOvXr+PSpUvIyMhAfn4+iouLodFomuT6hBBijp588kns3LkTs2fPhkwmw4ULFzBr1iz89NNPuHDhAi5evIjs7GxaX9SK0Zi1/7H2MWsAcOPGDTz77LMQi8U4ffp0ky8EXVZWhuzsbBQXF8POzg729vYQiUQQCoWwtbWFvb09xGIxxGIxJBIJhEIhjXMjVonGrFmPwsJCrFixAlu3boVOp4O3tzdGjhyJZ599FgEBAfD394dYLDZ1mKSZUbL2P5SsAYwx9OnTB4WFhRCJRAgMDMTrr7+OIUOGNOl1S0tLUVhYCJVKhfLycjDGwOfzuaRNn8CJxWLY29vD1taWS+LEYjG3rAshloqSNetz48YNLFy4EAkJCbC3t0fHjh0xevRoDBs2DIGBgdzqCMQ6CEwdADEfcXFxKCwsBACEhIRAIpHg3XffxfLly5s0YbO1tYWtrS33c3l5OUpLS7mHQqGARqMBn88Hn8/nkjWhUAiRSAQ7OzvY2dlxd+DEYjEEAvpoE0Jarg4dOmDdunU4ePAgFi9ejKSkJDx48ADHjx/HpEmT0L9/f7Rt25baOithFrcklEol3n//fXh5eUEikSAsLAxbt26t07EPHz7E5MmT4ebmBjs7O/Tp0wcHDx5s4ogt088//4ygoCAAQGZmJsLDwxEcHIwVK1Y067w/QqEQjo6OaNWqFXx9fREYGIj27dvD29sbjo6O0Gg0yM/PR1ZWFtLS0nDjxg1cvXoV165dw5UrV3D58mVcu3YNGRkZyMnJgUKhgFqtpskmSYNQ+0RMhcfjYfDgwdizZw/eeecdFBQUIDExEQsXLsTcuXPx3XffoU+fPrC3t0ePHj0QExNj6pBJEzGLlDwqKgqJiYlYvHgxOnTogM2bNyM6Oho6nQ4TJ06s9ji1Wo2IiAgUFhZi2bJl8PDwwA8//IBhw4YhPj4eAwcObMZX0fLduXMH06dPx+3bt/HgwQPMnz+fe65r167w9PREmzZtKj3atm3b5PMB8fl87g6aHmMMZWVlKCkpgUqlQlFRETQaDbRaLQQCAQQCAezs7LhuVJFIBKlUColEApFIBBsbG/D5/Cr/S2PjiB61T8TUxGIxXn/9dYwZMwbffPMN9u3bB7VajdOnT6NVq1aYOXMmTp48iYkTJ2LTpk2Iioqi4SEWxuRj1vbu3YuRI0dyDaDekCFDcOXKFaSnp1c77uLHH3/EW2+9hZMnT6JPnz4AAI1Gg65du0IqldZr6SQas1bxj5KzszNeeeUVHD16FOnp6UhISEBpaWmtd6UcHR0rJXD6/2/dunW1v8MDBw7g559/xp07d+Dv72+UMXKPdqOqVCqUlZVBq9WCz+fjypUriI+PR3Z2Nry8vDB69Gj06tULNjY2XJJmY2PDJXsCgYDrftXvU12CV1WyFxMTgwULFuDatWsICgrCzJkzERUV1ajXVxNLv15za4ntE41Zs3wXLlzAlClT4ODgAGdnZ7i6unJrkep0Ojz33HOQyWTccy4uLnBzc4ODgwOEQiH4fL5B+1aT5vgbt5RrNCWTJ2uvvvoqtm7dioKCAoO+9y1btmDixIk4ceIE+vbtW+WxkZGRyMjIwLVr1wy2L1q0CJ999hkyMzPh7e1dpzgoWatInN5991307dsX/fv3x7Fjx3Dy5EksX74c3bt3R3p6OjIzM5Geno6MjAzukZOTU+N5hUIhvLy8KiVxmZmZWLx4Mfr164ennnoKx48fx4kTJ5pkjJxWq0VsbCzmzJmDXr16ITQ0FFeuXMG1a9fw2muvoWvXrmCMgcfjVXrok7VHH/rt+n30yZo+URMKhUhISMDcuXPRo0cPPPnkkzh16hQSEhLw1VdfISIiotJ1AFR5/aq2V7Vt165dmDRpEiIiIjB48GDExcUhPj4eW7duxXPPPce9F4/fNXz050fPXZuYmBiMHTsWQ4YMwdChQ7F//34cOHAAO3bsaFGNYE1aYvtEyZp16NatG8LDw5GUlATGmEHyVd3fr0AggFQqhYODAxwcHODo6AiZTAZHR0c4OTnB2dmZS/BcXV1x7tw5vP766xg0aBAGDx6MAwcO4MCBA/j999+5xekbqznaEUtoq0yerPXp0wdarRZnzpwx2H7lyhV06dIFv/zyC1577bUqj/X09ET//v2xbds2g+179uzBqFGjsH///mr/0Ver1VCr1dzPCoUCbdq0sepkDahI2H755Rfcvn0bAQEBeP311xEZGVnjMSUlJcjMzOQSuUeTuros4q7/5gcA9+7dg1arha+vr9Fek15aWhr4fD68vb3B4/HAGKt0PcYYl3A9noQBqJSUVfX/erm5ubCxsYGLiwvXeObl5UGn08HDwwNA3ZIivdr2zc7Oho2NDTw8PLjX9/DhQ+h0Onh6etb5PI+qrnlgjCErKwt8Pp+73qeffoqff/4ZSqUSy5Ytq1f8dY3p0ffX29u7yacwaIntEyVr1kHfE/Ldd9/hxIkTyMnJwcaNG1FcXIzevXtDLpdDqVSiqKgIxcXFXJX94z0GteHxeHBycoJUKgWPx8ODBw+g0+ng5eUFoOo2ora04tHnMzMzuTWl9c/p55Sr65eZ2jx6DR6Ph1mzZuGHH36AQqHAypUrjXKNx/n4+Bi1+MPkY9by8vIQEBBQabv+H++8vLwaj62qfLkuxy5atAjz5s2rb7gAKj68AoEAcrm8Qcebs+7du+OXX34x2Jabm1vrcS4uLnBxcUFoaKjBdq1Wi5ycHGRlZSErKwv37t3j/qu/45CXl1fpd3Xjxo1GvpLqpaamVtpmzOs93j16584dgwQQALKysox+Z02/LS0trVLyk56ebrTX9yitVovMzEwAQEZGBgIDA7F79+4mu96jnJ2dmzxZa4ntE1BxB6WgoKDBxxPzN2nSJMycORNvvvkmevfujdOnTyMjIwMLFiyoNB6SMQaVSgW5XM49CgsLIZfLoVAoUFRUxCV1xcXFKCkpQWlpKYCK9iw/P7/S5yktLc1or0Wr1eLevXuVtle1rSF4PB50Oh2ysrIAVIzP7tChA3bt2oXbt28b5RqP8/DwsKxkDWjcN+6GHvvpp5/iww8/5H7Wf3OtC7FYjHbt2tFs0o0UHh4OR0dHzJ07F3K5HIwxfPHFF1Aqlfj666+Nfr1///vfcHR0xKxZs7g7T5ZyPcYY/vOf/0AqlWLmzJlc47RgwQIolUosXryY2+/R/z56/OP/X9NdNQCYO3cupFIpPvzwQ/B4PISEhGDDhg0QCoXw9/evV/UtY6zO+zPGIJFIIJVK63z+xmhp7ZOdnR0CAwOp+tnCderUCT4+Pvj222+xdu1adOjQARs3bsSoUaOMcn7GGAYOHAgHBwfMmzcPBQUF0Gq1WLx4MYqLi/Hll1/WeletLv+/YMECSKVSvPPOO9xz33//PYqLizFjxgyjvJbFixfD3t4eb7/9Nng8Hjp16oRNmzZBJBIZ9DoYi35Cd6NiJvbkk0+ynj17VtqenJzMALBffvml2mNbt27Nxo8fX2n77t27GQC2f//+Oschl8sZACaXy+t8DGmcHTt2MAAsMjKS/d///R+LjIxkAFhMTAxdj65nFqh9ItasOf7GLeUaTc3kydqrr77KpFIpKy8vN9i+ZcsWBoCdOHGi2mMjIyNZUFBQpe2LFi1iANi9e/fqHAc1hqaxY8cO1qNHD2ZnZ8d69OjR5H88dL2Wfb3mRu0TsXbN8TduKddoSiZP1vbu3csAsK1btxpsHzZsGPPy8mIajabaY3/88UcGgCUkJHDbysvLWefOnVnv3r3rFQc1hoSQx1H7RAgxByZP1hir+Abq7OzMVq5cyQ4dOsReffVVBoBt2rSJ22fq1KmMz+ezu3fvcttUKhXr3Lkza9OmDfv1119ZXFwcGzNmDBMIBOzIkSP1ioEaQ0JIVah9IoSYmlkUGMTExGDmzJmYPXs28vPzERQUhC1btmDChAncPlqtFlqt1mBgolgsxsGDB/Hxxx/jnXfeQUlJCcLCwrBv3z6aHZwQYhTUPhFCTM3k86yZi/pMOkkIIc2J2idCrBstHkYIIYQQYsYoWSOEEEIIMWNmMWbNHOh7gxUKhYkjIYQYm4ODQ72W2TI31D4RYrnq0j5RsvY/RUVFAFDnWcIJIS1HSx/rRe0TIZarLu0TFRj8j37dsJb+Ddyc6JfIycjIaNH/UJojem/rp6X/XdenfbLWzwa9but63YDlvHa6s1YPNjY28PHxMXUYFsnR0bFF/yGZM3pvrUND2idr/WzQ67Y+1vDaqcCAEEIIIcSMUbJGCCGEEGLGKFkjTUYsFmPOnDkQi8WmDsXi0HtLqmOtnw163db1ugHreu1UYEAIIYQQYsbozhohhBBCiBmjZI0QQgghxIxRskYIIYQQYsYoWSMNcuTIEfB4vCofCQkJBvueP38egwcPhlQqhZOTE6KionD79m0TRW5+ioqK8PHHH2PIkCFwd3cHj8fD3Llzq9y3Pu/lihUrEBQUBLFYDH9/f8ybNw/l5eVN+EqIKSmVSrz//vvw8vKCRCJBWFgYtm7dauqwjObQoUOYOnUqgoKCYG9vD29vb4wePRrnzp0z2G/y5MlVtktBQUEmirxxrLWtre73+Phrt7Tfd3VoUlzSKAsXLkR4eLjBti5dunD/f+3aNQwaNAhhYWHYtm0bVCoVZs+ejf79++PixYtwd3dv7pDNTl5eHlauXImuXbviueeew+rVq6vcrz7v5YIFCzBr1izMmDEDQ4YMQWJiIj7//HPcu3cPK1eubK6XRppRVFQUEhMTsXjxYnTo0AGbN29GdHQ0dDodJk6caOrwGu2nn35CXl4e3nvvPXTq1Ak5OTn4v//7Pzz55JPYv38/nn76aW5fW1tbHDp0yOB4W1vb5g7ZqKytrZ01axZef/31StufeeYZiMVi9OzZk9tmib/vShghDXD48GEGgG3fvr3G/caPH8/c3NyYXC7ntt29e5cJhUL28ccfN3WYLYJOp2M6nY4xxlhOTg4DwObMmVNpv7q+l7m5uUwikbDXXnvN4PgFCxYwHo/Hrly50jQvhJjMnj17GAC2efNmg+2RkZHMy8uLaTQaE0VmPA8ePKi0raioiLVq1YpFRERw215++WVmb2/fnKE1KWpr/3HkyBEGgH3++efcNkv7fVeHukFJk9FoNNi9ezfGjh1rsBSIr68vwsPDsXPnThNGZz70t+1rUp/3MjY2FiqVClOmTDE4x5QpU8AYwx9//GHU+Inp7dy5E1KpFOPHjzfYPmXKFGRlZeH06dMmisx4PDw8Km2TSqXo1KkTMjIyTBCR+bCWtnbNmjXg8XiYOnWqqUNpdpSskUZ56623IBAI4OjoiKFDh+L48ePcc7du3UJpaSlCQ0MrHRcaGoqbN29CpVI1Z7gtVn3ey+TkZABASEiIwX6enp5wc3PjnieWIzk5GcHBwRAIDEe26D8vlvo7l8vlOH/+PDp37mywvbS0FK1btwafz4ePjw/efvtt5OfnmyhK47D2tlYul+P3339HREQE/P39DZ6zxN/342jMGmkQmUyG9957D4MGDYKrqytu3ryJr7/+GoMGDcKePXswdOhQ5OXlAQBcXFwqHe/i4gLGGAoKCuDp6dnc4bc49Xkv8/LyIBaLYW9vX+W++nMRy5GXl4eAgIBK2/WfF0v9nb/11lsoLi7GzJkzuW1du3ZF165dufFcR48exbfffouDBw8iMTERUqnUVOE2CLW1FbZs2YLS0lJMmzbNYLul/b6rQ8kaaZBu3bqhW7du3M/9+/fHmDFjEBISgo8//hhDhw7lnqupi6+27j9iqK7vJb3n1sfafuezZs3Cr7/+ihUrVqBHjx7c9g8++MBgv8jISHTr1g3jxo3DqlWrKj1v7qitrbBmzRq4urpizJgxBtst7fddHeoGJUbj5OSEUaNG4dKlSygtLYWrqyuAqr/V5+fng8fjwcnJqZmjbJnq8166urpCpVKhpKSkyn2r+vZNWjZXV9dqPxtA1XdcWrJ58+Zh/vz5WLBgAd5+++1a9x8zZgzs7e0rTXXRUllbW3vp0iWcPXsWL774Yp3WAbW03zdAyRoxMva/pWZ5PB7atWsHW1tbXL58udJ+ly9fRmBgICQSSXOH2CLV573Uj1V7fN/s7Gzk5uYalPsTyxASEoKUlBRoNBqD7frPgCX9zufNm4e5c+di7ty5+Oyzz+p8HGMMNjaW80+eNbW1a9asAQC88sordT7G0n7flvNKiMkVFBRg9+7dCAsLg0QigUAgwDPPPIOYmBgUFRVx+6Wnp+Pw4cOIiooyYbQtS33ey2HDhkEikWD9+vUG51i/fj14PB6ee+65ZoqaNJcxY8ZAqVRix44dBts3bNgALy8v9O7d20SRGdeXX36JuXPn4vPPP8ecOXPqfNzvv/+OkpISPPnkk00YXfOxprZWrVZj06ZN6NWrV52/dFja7xsAzbNGGiY6Opp98sknbPv27ezw4cNs5cqVrGPHjkwgELC4uDhuv5SUFCaVStmAAQPY3r17WUxMDOvSpQvz8vJiDx8+NOErMC979+5l27dvZ2vXrmUA2Pjx49n27dvZ9u3bWXFxMWOsfu/l/PnzGY/HY5999hk7cuQI+/rrr5lYLGavvvqqKV4eaQaRkZHM2dmZrVy5kh06dIi9+uqrDADbtGmTqUMziqVLlzIAbNiwYezUqVOVHoxVzCvWt29ftnz5crZ37162b98+NmPGDCaRSFjnzp2ZUqk08auoP2tva7du3coAsJUrV1Z6zhJ/39WhZI00yKJFi1hYWBiTyWSMz+czd3d3NmbMGHbmzJlK+549e5ZFREQwOzs75ujoyJ577jl28+ZNE0Rtvnx9fRmAKh937tzh9qvPe7ls2TLWoUMHJhKJWNu2bdmcOXNYWVlZM70i0tyKiorYu+++y1q3bs1EIhELDQ1lW7ZsMXVYRjNw4MBq/0b09x3y8/PZmDFjmJ+fH7O1tWUikYi1b9+effzxx6ywsNDEr6BhrL2tjYyMZPb29kyhUFR6zhJ/39XhMfa/jm9CCCGEEGJ2aMwaIYQQQogZo2SNEEIIIcSMUbJGCCGEEGLGKFkjhBBCCDFjlKwRQgghhJgxStYIIYQQQswYJWuEEEIIIWaMkrVmoF/m5+zZs8163fz8fEyYMAEeHh60zND/HDlyBDweD0eOHOG27d27F3PnzjVZTI8aNGgQBg0axP189+5d8Hi8SktH1Wbz5s347rvvjBobIY1FbeE/SkpKMHfuXIO2qCGqatMmT54MPz+/ep0nKysLc+fOxcWLFxsVD2kaAlMHQJrOl19+iZ07d2Lt2rVo164dXFxcTB2SyXXv3h2nTp1Cp06duG179+7FDz/8YDYJ26M8PT1x6tQptGvXrl7Hbd68GcnJyXj//febJjBCWhBzbAtLSkowb948ADD4gmYMs2bNwnvvvVevY7KysjBv3jz4+fkhLCzMqPGQxqNkzYIlJyejXbt2mDRpUpNep7y8HDweDwKB+X+cHB0dm21x39LSUtja2jbqHGKx2LIWIybEBBrbFhrjb7k51ffLHTF/1A1qRlJTUzFx4kR4eHhALBYjODgYP/zwg8E+KpUKH330EcLCwiCTyeDi4oI+ffrgzz//5PbRd53Fx8cjJSUFPB6v0m3yx/n5+WHUqFHYuXMnQkNDIZFIEBAQgOXLlxvsp7/l/t///hcfffQRvL29IRaLcfPmTQDA2rVr0bVrV0gkEri4uGDMmDFISUmpdL1Vq1ahQ4cOEIvF6NSpEzZv3lzlrfuysjLMnz8fQUFBEIvFcHd3x5QpU5CTk1Nl/LGxsejevTtsbW0RFBSEtWvXVhm//r2YPHky9x7r3ycej4e7d+9y7/enn34Kf39/iEQieHt746233kJhYWGV14+JiUG3bt0gkUi4b81VYYzhq6++gq+vLyQSCbp37459+/ZV2q+qbtCcnBy89tpraNOmDfee9OvXD/Hx8QAqvqXv2bMHaWlpBq9Jb968eejduzdcXFzg6OiI7t27Y82aNXh85bm6vqcAcO/ePS4mkUgELy8vjBs3Dg8ePOD2USgU+Pe//23wXr7//vsoLi6u9n0i1skc2sKq/paTk5MxevRoODs7QyKRICwsDBs2bKh0jvT0dLz44osG8f/f//0fdDodF5e7uzuAir9HfVyTJ0+u8X25du0ahg0bBjs7O7i5ueH1119HUVFRpf2qaku3b9+O3r17QyaTwc7ODgEBAZg6dSqAinaxZ8+eAIApU6Zw8eh7G86ePYsJEybAz88Ptra28PPzQ3R0NNLS0gyuoe/mPnz4MN544w24ubnB1dUVUVFRyMrKqhTn5s2b0adPH0ilUkilUoSFhWHNmjUG+8THxyMiIgKOjo6ws7NDv379cPDgwRrfJ4tk2qVJrcO6desYAJaYmFjtPleuXGEymYyFhISwjRs3sgMHDrCPPvqI2djYsLlz53L7FRYWssmTJ7P//ve/7NChQyw2Npb9+9//ZjY2NmzDhg2MMcZUKhU7deoU69atGwsICGCnTp1ip06dYnK5vNrr+/r6Mm9vb9a2bVu2du1atnfvXjZp0iQGgH399dfcfocPH2YAmLe3Nxs3bhz766+/2O7du1leXh5buHAhA8Cio6PZnj172MaNG1lAQACTyWTsxo0b3Dl++eUXBoCNHTuW7d69m/3666+sQ4cOzNfXl/n6+nL7abVaNmzYMGZvb8/mzZvH4uLi2OrVq5m3tzfr1KkTKykpMYjfx8eHderUiW3cuJHt37+fjR8/ngFgR48erRT/4cOHGWOM3bx5k40bN44B4N6nU6dOMZVKxXQ6HRs6dCgTCARs1qxZ7MCBA2zp0qXM3t6edevWjalUKoPre3p6soCAALZ27Vp2+PDhKhda1pszZw4DwKZNm8b27dvHVq5cyby9vVnr1q3ZwIEDuf3u3LnDALB169Zx24YOHcrc3d3ZypUr2ZEjR9gff/zBZs+ezbZu3cp9lvr168dat25t8Jr0Jk+ezNasWcPi4uJYXFwc+/LLL5mtrS2bN29epc9EXd7TzMxM5unpydzc3Ng333zD4uPj2W+//camTp3KUlJSGGOMFRcXs7CwMIN9li1bxmQyGXv66aeZTqer9r0ilqOltIVV/S1fu3aNOTg4sHbt2rGNGzeyPXv2sOjoaAaALVmyhDv+4cOHzNvbm7m7u7Off/6ZxcbGsrfffpsBYG+88QYXV2xsLNcG6OOqadH17Oxs5uHhwby9vdm6deu4Nrpt27YGbRpjjL388ssGbenJkycZj8djEyZMYHv37mWHDh1i69atYy+99BJjjDG5XM79bj7//HMunoyMDMYYY9u3b2ezZ89mO3fuZEePHmVbt25lAwcOZO7u7iwnJ6fS7zcgIIC98847bP/+/Wz16tXM2dmZhYeHG7yeWbNmMQAsKiqKbd++nR04cIB98803bNasWdw+//3vfxmPx2PPPfcci4mJYbt27WKjRo1ifD6fxcfHV/teWSJK1ppBXRqooUOHMh8fn0qNyNtvv80kEgnLz8+v8jiNRsPKy8vZtGnTWLdu3QyeGzhwIOvcuXOdYvT19WU8Ho9dvHjRYHtkZCRzdHRkxcXFjLF/kp0BAwYY7FdQUMBsbW3ZiBEjDLanp6czsVjMJk6cyBirSMBat27NevfubbBfWloaEwqFBg3Mli1bGAC2Y8cOg30TExMZAPbjjz8axC+RSFhaWhq3rbS0lLm4uLDp06dz2x5P1hhj7K233mJVfW/RN6ZfffWVwfbffvuNAWArV640uD6fz2fXr1+vdJ7HFRQUMIlEwsaMGWOw/cSJEwxArcmaVCpl77//fo3XGDlypMF7WR2tVsvKy8vZF198wVxdXQ2Sprq+p1OnTmVCoZBdvXq12ussWrSI2djYVPob+P333xkAtnfv3lpjJS1fS2kLq/pbnjBhAhOLxSw9Pd1g+/Dhw5mdnR0rLCxkjDE2Y8YMBoCdPn3aYL833niD8Xg87rw5OTkMAJszZ06d4vrkk0+qbaNrS9aWLl3KAHAxVkXfrj7a1lRHo9EwpVLJ7O3t2bJly7jt+t/vm2++abD/V199xQCw+/fvM8YYu337NuPz+WzSpEnVXqO4uJi5uLiwZ555xmC7VqtlXbt2Zb169ao1TktC3aBmQKVS4eDBgxgzZgzs7Oyg0Wi4x4gRI6BSqZCQkMDtv337dvTr1w9SqRQCgQBCoRBr1qypsrvxUTqdzuDcWq3W4PnOnTuja9euBtsmTpwIhUKB8+fPG2wfO3aswc+nTp1CaWlppdv4bdq0wdNPP83dtr5+/Tqys7Px/PPPG+zXtm1b9OvXz2Db7t274eTkhGeeecYg7rCwMLRu3bpSV0ZYWBjatm3L/SyRSNChQ4dKt+rr6tChQwBQ6TWNHz8e9vb2lW7Fh4aGokOHDrWe99SpU1CpVJXGz/Tt2xe+vr61Ht+rVy+sX78e8+fPR0JCAsrLy2s95lGHDh3C4MGDIZPJwOfzIRQKMXv2bOTl5eHhw4cG+9blPd23bx/Cw8MRHBxc7TV3796NLl26ICwszOB3OXTo0Fq7pYj1MJe2sKq/5UOHDiEiIgJt2rQx2D558mSUlJTg1KlT3H6dOnVCr169Ku3HGOPalfo6fPhwtW10bfRdnM8//zy2bduGe/fu1evaSqUSn3zyCQIDAyEQCCAQCCCVSlFcXFzle/3ss88a/BwaGgoAXLsRFxcHrVaLt956q9prnjx5Evn5+Xj55ZcNflc6nQ7Dhg1DYmKiVQ2hoGTNDOTl5UGj0WDFihUQCoUGjxEjRgAAcnNzAQAxMTF4/vnn4e3tjU2bNuHUqVNITEzE1KlToVKparzO1KlTDc4dERFh8Hzr1q0rHaPflpeXZ7Dd09Oz0muoajsAeHl5cc/r/9uqVatK+z2+7cGDBygsLIRIJKr0vmRnZ3PviZ6rq2ulc4rFYpSWllbaXhd5eXkQCATc2BI9Ho+H1q1b1/qe1HReoOb3uya//fYbXn75ZaxevRp9+vSBi4sL/vWvfyE7O7vWY8+cOYMhQ4YAqBg3eOLECSQmJmLmzJkAUOm9qst7mpOTAx8fnxqv++DBA1y6dKnS79HBwQGMsUq/S2KdzKUtrOpvOS8vr9r2Tf98ffarr7y8vAa3GQMGDMAff/wBjUaDf/3rX/Dx8UGXLl2wZcuWOl174sSJ+P777/HKK69g//79OHPmDBITE+Hu7l5l+/p4uyEWiwH8077oxxzX1G7ox7uOGzeu0mdhyZIlYIwhPz+/TvFbAvMv37MCzs7O4PP5eOmll6r9puHv7w8A2LRpE/z9/fHbb78ZDBpXq9W1Xmfu3Ll4++23uZ8dHBwMnq/qH3v9tsf/+B699qPP379/v9I5srKy4ObmZrDfowPPq7u+fnBqbGxsla/n8fiNzdXVFRqNBjk5OQYJG2MM2dnZ3LdVvcffk5rOC1T/ftc2P5Kbmxu+++47fPfdd0hPT8dff/2FGTNm4OHDh9W+V3pbt26FUCjE7t27IZFIuO1//PFHnWKviru7OzIzM2uN2dbWtsriBP3zhJhLW1jV37Krq2u17RsAgzauLvvVl6ura41tdG1Gjx6N0aNHQ61WIyEhAYsWLcLEiRPh5+eHPn36VHucXC7H7t27MWfOHMyYMYPbrlarG5ws6dvTzMzMSncq9fTv04oVK6qtiK/qS7+lomTNDNjZ2SE8PBwXLlxAaGgoRCJRtfvyeDyIRCKDxiQ7O9ugAqo6fn5+NSYCV65cQVJSksFt9s2bN8PBwQHdu3ev8dx9+vSBra0tNm3ahPHjx3PbMzMzcejQIYwbNw4A0LFjR7Ru3Rrbtm3Dhx9+yO2Xnp6OkydPct8+AWDUqFHYunUrtFotevfuXevra6hHv/U9Wp4fERGBr776Cps2bcIHH3zAbd+xYweKi4srfRuvqyeffBISiQS//vqrQXfyyZMnkZaWVq/JLNu2bYu3334bBw8exIkTJwxeU1XfePVTrPD5fG5baWkp/vvf/zbotQDA8OHD8d///hfXr19Hx44dq9xn1KhRWLhwIVxdXbl/bAl5nLm0hVWJiIjAzp07kZWVZdBObdy4EXZ2dlxCERERgUWLFuH8+fMG7ebGjRvB4/EQHh4OoPLdptqEh4fjq6++qrKNrg+xWIyBAwfCyckJ+/fvx4ULF9CnT59q4+HxeGCMcc/rrV69ulL3cV0NGTIEfD4fP/30U7WJYr9+/eDk5ISrV68aJNbWipK1ZnTo0CFuSohHjRgxAsuWLcNTTz2F/v3744033oCfnx+Kiopw8+ZN7Nq1ixvnoC8pf/PNNzFu3DhkZGTgyy+/hKenJ1JTUxsVn5eXF5599lnMnTsXnp6e2LRpE+Li4rBkyRLY2dnVeKyTkxNmzZqFzz77DP/6178QHR2NvLw8zJs3DxKJBHPmzAEA2NjYYN68eZg+fTrGjRuHqVOnorCwEPPmzYOnpydsbP7pmZ8wYQJ+/fVXjBgxAu+99x569eoFoVCIzMxMHD58GKNHj8aYMWMa9ZoBICQkBACwZMkSDB8+HHw+H6GhoYiMjMTQoUPxySefQKFQoF+/frh06RLmzJmDbt264aWXXmrQ9ZydnfHvf/8b8+fPxyuvvILx48cjIyMDc+fOrbVLQy6XIzw8HBMnTkRQUBAcHByQmJiI2NhYREVFGbymmJgY/PTTT+jRowdsbGzwxBNPYOTIkfjmm28wceJEvPbaa8jLy8PSpUsrNcT18cUXX2Dfvn0YMGAAPvvsM4SEhKCwsBCxsbH48MMPERQUhPfffx87duzAgAED8MEHHyA0NBQ6nQ7p6ek4cOAAPvrooyZNyIl5Mfe2sCpz5szB7t27ER4ejtmzZ8PFxQW//vor9uzZg6+++goymQwA8MEHH2Djxo0YOXIkvvjiC/j6+mLPnj348ccf8cYbb3Bj4RwcHODr64s///wTERERcHFxgZubW7VJ5Pvvv4+1a9di5MiRmD9/Plq1aoVff/0V165dqzX22bNnIzMzExEREfDx8UFhYSGWLVsGoVCIgQMHAqiYm83W1ha//vorgoODIZVK4eXlBS8vLwwYMABff/01F9/Ro0exZs0aODk5Nei99PPzw2effYYvv/wSpaWliI6Ohkwmw9WrV5Gbm4t58+ZBKpVixYoVePnll5Gfn49x48bBw8MDOTk5SEpKQk5ODn766acGXb9FMml5g5XQV8hU97hz5w5jrKLyb+rUqczb25sJhULm7u7O+vbty+bPn29wvsWLFzM/Pz8mFotZcHAwW7VqFTcVxKPqWwE1cuRI9vvvv7POnTszkUjE/Pz82DfffGOwn76acvv27VWeZ/Xq1Sw0NJSJRCImk8nY6NGj2ZUrVyrtt3LlShYYGMhEIhHr0KEDW7t2LRs9enSlKq7y8nK2dOlS1rVrVyaRSJhUKmVBQUFs+vTpLDU1tVL8jxs4cKBBdWVV1aBqtZq98sorzN3dnfF4PIPfSWlpKfvkk0+Yr68vEwqFzNPTk73xxhusoKCgyvevrnQ6HVu0aBFr06YNE4lELDQ0lO3atatSvI9Xg6pUKvb666+z0NBQ5ujoyGxtbVnHjh3ZnDlzuIpdxhjLz89n48aNY05OTtxr0lu7di3r2LEjE4vFLCAggC1atIitWbPG4HXX5z1ljLGMjAw2depU1rp1ayYUCpmXlxd7/vnn2YMHD7h9lEol+/zzz1nHjh25z0dISAj74IMPWHZ2dp3fO9JytaS2sCqXL19mzzzzDJPJZEwkErGuXbtWWT2ZlpbGJk6cyFxdXZlQKGQdO3ZkX3/9NdNqtQb7xcfHs27dujGxWMwAsJdffrnG2K5evcoiIyOZRCJhLi4ubNq0aezPP/+stRp09+7dbPjw4czb25uJRCLm4eHBRowYwY4dO2Zw/i1btrCgoCAmFAoNKlUzMzPZ2LFjmbOzM3NwcGDDhg1jycnJzNfX1yDm6qp9q2p3GWNs48aNrGfPnlzb3q1bt0rv59GjR9nIkSOZi4sLEwqFzNvbm40cObLaf4MsFY+xx2bCJFbJz88PXbp0we7du01y/cLCQnTo0AHPPfccVq5caZIYCCGEEHNE3aCk2WVnZ2PBggUIDw+Hq6sr0tLS8O2336KoqKje69kRQgghlo6SNdLsxGIx7t69izfffBP5+fnc4Nyff/4ZnTt3NnV4hBBCiFmhblBCCCGEEDNGk+ISQgghhJgxStYIIYQQQswYJWuEEEIIIWaMkjVCCCGEEDNGyRohhBBCiBmjZI0QQgghxIxRskYIIYQQYsYoWSOEEEIIMWOUrBFCCCGEmLH/B5HYUeKbEmASAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, 2, figsize=(7, 3))\n", "\n", "max_val = 130\n", "n_hist = 8\n", "xrange = np.arange(0, max_val, int(max_val/n_hist))+int(max_val/n_hist)/2\n", "dist = np.histogram(branch_len+1, xrange-int(max_val/n_hist)/2)\n", "th_dist = [spx.est_mr.BranchLength(*az.summary(idata_bl)['mean']).prob(i) for i in xrange]\n", " \n", "ax[0].plot(xrange[:-1], dist[0]/sum(dist[0]), 'o-', mfc='none', c='black', ms=4, label='Data')\n", "ax[0].plot(xrange, th_dist/sum(th_dist), c='tab:grey',alpha=0.3, label='Estimated')\n", "ax[0].fill_between(xrange, th_dist/sum(th_dist), color='tab:grey', alpha=0.3)\n", "ax[0].spines['right'].set_visible(False)\n", "ax[0].spines['top'].set_visible(False)\n", "ax[0].set_xlabel('Leaf-progenitor distance')\n", "ax[0].set_ylabel('Density')\n", "ax[0].legend()\n", "print(f'lp_{chisquare(f_obs = dist[0]/sum(dist[0]), f_exp=th_dist[:-1]/sum(th_dist[:-1]), ddof=2)}')\n", "\n", "\n", "\n", "alpha = n_stem[:,-1]+n_nonstem[:,-1]\n", "alpha = alpha / sum(alpha)\n", "\n", "max_val = 100\n", "n_hist = 10\n", "mutdist = np.histogram(mutnum, np.arange(0, max_val, int(max_val/n_hist)))\n", "th_dist = np.zeros(max_val)\n", "for i, a in enumerate(alpha):\n", " th_dist = th_dist + a*poisson((i+1)*2).pmf(range(max_val))\n", "th_dist_x = []\n", "th_dist_y = []\n", "for i in range(n_hist):\n", " th_dist_x.append(i*int(max_val/n_hist))\n", " th_dist_y.append(np.sum(th_dist[i*int(max_val/n_hist):(i+1)*int(max_val/n_hist)]))\n", "ax[1].plot(mutdist[1][:-1], mutdist[0]/sum(mutdist[0]), 'o-', mfc='none', c='black', ms=4, label='Data')\n", "ax[1].plot(th_dist_x, th_dist_y, c='tab:grey',alpha=0.3, label='Estimated')\n", "ax[1].fill_between(th_dist_x, th_dist_y, color='tab:grey', alpha=0.3)\n", "ax[1].spines['right'].set_visible(False)\n", "ax[1].spines['top'].set_visible(False)\n", "ax[1].set_xlabel('Leaf-root distance')\n", "ax[1].set_ylabel('Density')\n", "ax[1].legend()\n", "print(f'lr_{chisquare(f_obs = mutdist[0]/sum(mutdist[0]), f_exp=th_dist_y[:-1]/sum(th_dist_y[:-1]), ddof=5)}')" ] }, { "cell_type": "code", "execution_count": null, "id": "1561fd2b-fa45-41ae-9278-5a31e850eb96", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "43ba76a0-d0bf-4dd2-b2df-57b51c453af2", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "82606b06-9126-4466-86e5-4ab46d980e2d", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.12" } }, "nbformat": 4, "nbformat_minor": 5 }